\documentclass[12pt]{article}
\input{../Bibliography/commoncommands}

\usepackage{tikz,tikz-3dplot}
\usetikzlibrary{arrows}

\newenvironment{myfont}{\fontfamily{\ttdefault}\selectfont}{\par}

 \mathchardef\mhyphen="2D

\begin{document}
\bibliographystyle{unsrt}

\vspace{1.0in}

% To compile geom_notes:
%
% 1. Make sure you updated your repo and compiled the latest produciton version of fds.
%
% 2. Go to fds/Manuals/FDS_User_Guide/
%
% 3. Type:
%
% [mnv@burn FDS_User_Guide]$ ./make_geom_figures.sh -J
%
% The -J is for the case where you compiled the target in fds/Build/impi_intel_linux_64 (ifort + intel mpi). If you don't use this flag,
% it assumes running with the target in fds/Build/mpi_intel_linux_64 (ifort + openmpi).
% This script runs all the cases necessary to make smokeview figures, runs smokeview in batch mode and makes these figs. They are saved in
% fds/Manuals/FDS_User_Guide/SCRIPT_FIGURES.
%
% 4. Use your latex editor to compile geom_notes.tex or type:
%
% [mnv@burn FDS_User_Guide]$ ./make_geom_notes.sh
%
% You should get the geom_notes.pdf in the same directory.
%\chapter{Building the Model (FIXME UG Chapter 7)}
%
%\textbf{FIXME}: [EG] in the following paragraph modifications are \underline{underlined}
%
%A considerable amount of work in setting up a calculation lies in specifying the geometry of the space to
%be modeled and applying boundary conditions to the solid surfaces. The geometry is described in terms of
%\underline{solids} that can heat up, burn, conduct heat, etc.; and vents from which air or fuel can be
%either injected into, or drawn from, the flow domain. A boundary condition needs to be assigned to each
%\underline{solid surface} and vent describing its thermal properties. A fire is just one type of boundary condition.
%This chapter describes how to build the model.
%
%\textbf{FIXME}: the following text is written following the structure of the {\ct OBST} namelist section

\section{Creating Unstructured Solids: The \texorpdfstring{{\ct GEOM}}{GEOM} Namelist Group (Table \ref{tbl:GEOM})}
\label{info:GEOM}

The namelist {\ct GEOM}\ describes one or more unstructured closed geometric surfaces that enclose solid portions
of the fluid domain. These surfaces consist of a collection of triangular faces, where each face is built
from three vertices. The user can assign a specific boundary condition to each of the faces.

\subsection{Basics}
\label{info:GEOM_Basics}

A simple form of the {\ct GEOM}\ namelist defining an unstructured solid is given by

\begin{verbatim}
&GEOM ID='UNSTRUCTURED_SOLID'
      SURF_ID='FIRE','INERT'
      VERTS= -1.0, -1.0,  0.0,
              1.0, -1.0,  0.0,
              0.0,  1.0,  0.0,
              0.0,  0.0,  1.0,
      FACES= 1,3,2, 2,
             1,4,3, 1,
             3,4,2, 1,
             2,4,1, 0,
/
\end{verbatim}

\noindent where {\ct ID}\ specifies the solid name, in this case {\ct UNSTRUCTURED\_SOLID},
{\ct VERTS}\ specifies the $(x_1,y_1,z_1, \quad x_2,y_2,z_2, \quad ...)$ coordinates of the vertices,
and {\ct FACES}\ specifies a list of triangular faces and the {\ct SURF\_ID}\ that should be
assigned to each of them.

A triangular face is defined by 3 vertex indices and one boundary condition index: $(v_1,v_2,v_3,b)$.
Each vertex index $(v_i)$ ranges from 1 to the number of vertices found in the {\ct VERTS} parameter.
The order of the 3 vertex indices in the face defines which side of the triangle is considered facing
the outside of the solid. Observed from the outside, the vertices of each triangle should always appear
to have a counter-clockwise winding order, so that their order rotates counter-clockwise around the
triangle's center and the face normal vector is oriented according to the right-hand rule.

The boundary condition index $b$ ranges from 0 to the number of boundary condition {\ct IDs} found in
the {\ct SURF\_ID} parameter. A value greater than 0 designates, through the {\ct SURF\_ID} list, which
{\ct SURF} line (Section~\texttt{info:SURF}) to apply at the selected face. This index $b$ is always defined locally to the
{\ct GEOM SURF\_ID} surface list. A value equal to 0 applies the \textit{default} boundary condition to the selected face.
Therefore, in the above example, the first face has an index $b=2$ which corresponds to the {\ct 'INERT'} surf ID, the
second and third faces have $b=1$ and are of type {\ct 'FIRE'}, and the fourth will be of \textit{default} surface ID.

\subsection{Triangulated surfaces quality} \label{triangulated_surfaces_quality}

In order for FDS to correctly detect the solid portion of the volume from the rest
of the fluid domain, some tests are performed at start on the quality of
the triangulated surfaces.

The following conditions are enforced by FDS on the collection of vertices and faces of each {\ct GEOM}\ namelist:

\begin{enumerate}
\item Each vertex must be unique, shared by edges, be part of an only well-formed surface so that its neighbor ring of faces could in theory be continuously deformed to a disk: a vertex that abides to these conditions is called \textit{non-degenerate} and \textit{manifold}. Examples of good and bad vertices are presented in Figure \ref{figure:GEOM_non_manifold_vert}.

\item Each edge connecting vertices must have a non-zero length, and must always be shared by exactly two faces. An edge that abides to these conditions is called \textit{non-degenerate} and \textit{manifold}. See examples of good and bad edges in Figure \ref{figure:GEOM_non_manifold_edge}.

\item Each face cannot intersect the others and must have a non-zero area; such a face is called \textit{non-intersecting} and \textit{non-degenerate}.

\item The vertex orderings for the faces shall be chosen so that adjacent faces
have outward consistent normals, as illustrated in Figure \ref{figure:GEOM_oriented_surface}.

\item Each {\ct GEOM}\ namelist may define one or more distinct volumes (\textit{unconnected} volumes) that should never mutually intersect or self-intersect. See an example of a bad intersecting geometry in Figure \ref{figure:GEOM_self_intersecting}.

\item An upper limit is set to the number of triangles for each {\ct MESH} cell; this condition reasonably limits the complexity of the discretization that can currently be elaborated by the immersed boundary method solver.

\end{enumerate}

If the previous conditions are respected, a well-formed {\ct GEOM}\ namelist is obtained that represents one or more \textit{closed}, \textit{manifold}, \textit{orientable}, \textit{non intersecting} triangulated surfaces, that enclose well-defined volumes. An example of a well-formed triangulated surface is presented in Figure \ref{figure:GEOM_ok_surface_2}.

\pagebreak
Whenever one of the checks on these conditions fails, FDS outputs an error message describing the problem and its location in the domain. For example:\nopagebreak
\begin{verbatim}
ERROR: GEOM ID='Cube': Non manifold geometry at edge: ...
\end{verbatim}

\begin{figure}
	\centering
	\begin{minipage}{.45\textwidth}
		\centering
		\includegraphics[width=1\textwidth]{../../../fig/fds/Geometry_Figures/GEOM_oriented_surface.png}
		\caption{\textbf{Good {\ct GEOM}:} consistent normals on a manifold orientable surface}
		\label{figure:GEOM_oriented_surface}
	\end{minipage}%
	\hfill
	\begin{minipage}{.45\textwidth}
		\centering
		\includegraphics[width=1\textwidth]{../../../fig/fds/Geometry_Figures/GEOM_ok_surface_2.png}
		\caption{\textbf{Good {\ct GEOM}:} a well formed triangulated surface}
		\label{figure:GEOM_ok_surface_2}
	\end{minipage}%
\end{figure}

\begin{figure}
	\centering
	\begin{minipage}{.45\textwidth}
		\centering
		\includegraphics[width=1\textwidth]{../../../fig/fds/Geometry_Figures/GEOM_non_manifold_vert.png}
		\caption{\textbf{Bad {\ct GEOM}:} the vertex shared by the two cubes is not manifold, because it is connected to two rings of faces, that cannot be continuously deformed to a single disk. The other vertices are manifold.}
		\label{figure:GEOM_non_manifold_vert}
	\end{minipage}%
	\hfill
	\begin{minipage}{.45\textwidth}
		\centering
		\includegraphics[width=1\textwidth]{../../../fig/fds/Geometry_Figures/GEOM_non_manifold_edge.png}
		\caption{\textbf{Bad {\ct GEOM}:} the edge shared by the two cubes is not manifold, because more than two faces are incident on it. The other edges are manifold.\newline}
		\label{figure:GEOM_non_manifold_edge}
	\end{minipage}%
\end{figure}

\begin{figure}
	\centering
	\begin{minipage}{.45\textwidth}
		\centering
		\includegraphics[width=1\textwidth]{../../../fig/fds/Geometry_Figures/GEOM_open_surface.png}
		\caption{\textbf{Bad {\ct GEOM}:} FDS cannot differentiate the internal volume from the rest of the domain, because the triangulated surface is open}
		\label{figure:GEOM_open_surface}
	\end{minipage}%
	\hfill
	\begin{minipage}{.45\textwidth}
		\centering
		\includegraphics[width=1\textwidth]{../../../fig/fds/Geometry_Figures/GEOM_open_thick_surface.png}
		\caption{\textbf{Good {\ct GEOM}:} by adding a thickness FDS can differentiate the internal volume from the rest of the domain\newline}
		\label{figure:GEOM_open_thick_surface}
\end{minipage}%
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=.5\textwidth]{../../../fig/fds/Geometry_Figures/GEOM_self_intersecting.png}
	\caption{\textbf{Bad {\ct GEOM}:} This volume is self-intersecting}
	\label{figure:GEOM_self_intersecting}
\end{figure}

\subsection{Intersections and interactions}

Several types of intersections can happen between {\ct GEOM}s and between {\ct GEOM}s and other namelists:

\begin{itemize}
	\item {\ct GEOM} - {\ct GEOM} and {\ct GEOM} - {\ct OBST}: whenever these types of intersection happen, FDS performs a geometric boolean operation uniting the two volumes, respecting the boundary conditions of the surfaces. If the geometry is too complex, this operation can fail. In that case FDS will print out a warning to the screen (or to standard error) and fill the corresponding cell with a solid condition;
	\item {\ct GEOM} - {\ct HOLE}: currently there is no interaction.
	\item {\ct GEOM} - {\ct VENT}: it is currently not possible to set boundary condition patches on {\ct GEOM} surfaces with {\ct VENT} namelists.
\end{itemize}

Currently no control logic is available on {\ct GEOM} namelists.

\subsection{Coloring}

FIXME Colors and texture maps

\subsection{Self-generated geometries}

The {\ct GEOM}\ namelist allows the quick definition of particular kinds of geometric objects, blocks using {\ct XB}, spheres using {\ct SPHERE\_ORIGIN}\ and {\ct SPHERE\_RADIUS}
and a 2D terrain elevations using {\ct ZVALS}. FDS generates vertices and faces to represent these objects, equivalent to what would have been defined using {\ct VERTS}\ and {\ct FACES}.

\subsubsection{Blocks}
A {\ct GEOM}\ namelist defining a block is given by

\begin{verbatim}
&GEOM ID='block'
      SURF_ID='S1'
      XB=0.0,1.0,0.0,1.0,0.0,1.0
/
\end{verbatim}

\noindent where {\ct XB=xmin,xmax,ymin,ymax,zmin,zmax} defines the min and max bounds of the block.
The {\ct XB}\ parameter is used in the same way as on an {\ct OBST}\ or {\ct VENT}\ line.
A block may be refined into many parts by specifying the {\ct IJK} parameter.
For example, {\ct IJK=8,6,4} would split the block into 8 parts along the $x$ dimension,
6 parts along the $y$ dimension and 4 parts along the $z$ dimension.
By default, blocks are discretized so that the block faces are consistent in size with the grid resolution.

The {\ct SURF\_ID}\ parameter assigns the specified boundary condition to all the generated faces.

\subsubsection{Spheres}
A {\ct GEOM}\ namelist defining a sphere centered at $(0,0,0)$ with radius $1$ is given by

\begin{verbatim}
&GEOM ID='sphere'
      SURF_ID='S1'
      SPHERE_ORIGIN=0.0,0.0,0.0, SPHERE_RADIUS=1.0
/
\end{verbatim}

\noindent Spheres are discretized by default so that each face is consistent in size with the grid resolution.
One may specify a {\ct N\_LEVEL}\ parameter which defines the number of times the sphere is split.


\subsubsection{Cylinders}
\label{info:GEOM_cyls}

A {\ct GEOM} namelist defining a cylinder geometry with centroid at {\ct CYLINDER\_ORIGIN}$=(0,0,0)$, axis {\ct CYLINDER\_AXIS}$=(1.,0,0)$, length {\ct CYLINDER\_LENGTH}$=2.$ and radius {\ct CYLINDER\_RADIUS}$=0.5$ is given by

\begin{verbatim}
&GEOM ID='cylinder', SURF_ID='S1',
              CYLINDER_LENGTH=2.,
              CYLINDER_RADIUS=.5,
              CYLINDER_ORIGIN=0.,0.,0.,
              CYLINDER_AXIS=1.,0.,0.,
              CYLINDER_NSEG_AXIS=6,
              CYLINDER_NSEG_THETA=48 /
 \end{verbatim}

where the cylinder is discretized using 6 segments along its axis and 48 segments across its circumference. The keyword {\ct SURF\_ID} can be replaced by {\ct SURF\_IDS='S1','S2','S3'}, where the strings {\ct 'S1','S2','S3'} correspond to the surface IDs of top, sides and bottom surfaces respectively.


\subsubsection{2D Terrain elevations}
A {\ct GEOM}\ namelist defining a 2D terrain elevation is given by

\begin{verbatim}
&GEOM ID='terrain'
      SURF_ID='S1'
      IJK=ivals,jvals, XB=xmin,xmax,ymin,ymax,
      EXTEND_TERRAIN=etlog, ZVAL_HORIZON=zvhr,
      ZVALS=...
/
\end{verbatim}

\noindent where {\ct XB}\ defines a rectangular region bounded by $(x_{\rm min}, y_{\rm min})$  and $(x_{\rm max}, y_{\rm max})$
and where {\ct ZVALS} defines elevation data in this region.
{\ct IJK}\ specifies how many vertices occur in this region along each dimension.
In this example {\ct ivals}\ values occur along the x dimension and {\ct jval}\ values occur along the y dimension.
The {\ct ZVALS} keyword is used to specify elevations at each $(x,y)$ location.
The elevation data specified after the {\ct ZVALS}\ keyword is arranged in row major order.
The first row contains {\ct ivals} elevation values occurring at the the ymax position from xmin to xmax.
Similarly, the last row contains
{\ct ivals}\ elevation values for the ymin position again from xmin to xmax.
There are then {\ct jvals}\ rows and {\ct ivals}\ columns of elevation data.
The optional fields {\ct EXTEND\_TERRAIN} and {\ct ZVAL\_HORIZON} are used in the following form. If the computational domain is set larger than the {\ct XB} defined terrain patch, such that it contains the terrain completely, setting {\ct etlog=.TRUE.} allows for the terrain to be extended to the domain boundaries. This can be helpful when trying to reduce boundary effects on a calculation.
If, additionally, the real number {\ct zvhr} is defined for {\ct ZVAL\_HORIZON} on the {\ct \&GEOM} line, this elevation will be used on the domain boundaries when extending the domain.
As with the blocks and spheres, FDS uses the information provided by these keywords to construct
vertices and triangular faces.

Additionally, terrain geometries can be specified defining a triangulation with a convex regular boundary, as seen from the top. This triangulation will represent the terrain of interest and will define the top boundary of a terrain geom. The terrain surface is specified by {\ct VERTS} and {\ct FACES} as in a basic geometry (section~\ref{info:GEOM_Basics}), and the field {\ct IS\_TERRAIN} needs to be set to {\ct TRUE}, such that FDS knows this is a terrain unstructured surface. Surface types per faces are specified in the same manner as in basic geometries. Also, the {\ct EXTEND\_TERRAIN} logic can also be applied.

\subsubsection{Geometries extruded from simple planar polygons}

In this case, a geometry can be defined specifying the {\ct VERTS}, connectivity {\ct POLY}, and extrusion distance {\ct EXTRUDE}. For example:

\begin{verbatim}
&GEOM ID='ExtPoly', SURF_ID='Poly1',
VERTS =
   -1.0, -1.0,  1.0,
    0.5, -1.0,  1.75,
    1.0, -1.0,  2.0,
    0.0,  0.0,  1.5,
    1.0,  1.0,  2.0,
   -1.0,  1.0,  1.0,
POLY = 1, 2, 4, 3, 5, 6, 
EXTRUDE=0.5 /
\end{verbatim}

The polygon defined by {\ct POLY} must be simple, non degenerate (no segment intersections, no repeated vertices) and all segments must lay on the same plane in three dimensional space. It can have a maximum of 1000 vertices. The value of {\ct EXTRUDE} can be a positive or negative real, but not zero. A single surface ID can be provided, or top, sides and bottom surface IDs through a  {\ct SURF\_IDS} field, in the same manner as the cylinders of section~\ref{info:GEOM_cyls}. Note that the bottom surface in {\ct SURF\_IDS} always refers to the polygon defined in {\ct POLY}, from where the extrusion operation is performed.


\subsection{Generating complex geometries}

Generating well-formed complex geometries is a time-consuming and error-prone process.
Several third-party pre-processing tools are provided by the FDS-SMV community for automatically translating CAD files to {\ct GEOM}\ namelists.

These tools should always enforce quality checks on the triangulated surfaces and produce well-formed {\ct GEOM}\ namelists.
Links to these tools can be found at the project home page at https://pages.nist.gov/fds-smv/.



%\subsection{Defining Surfaces}
%
%\subsubsection{General Surfaces}
%
%Unstructured geometric surfaces are defined using the {\ct GEOM}\ namelist.
%These surfaces consist of a collection of triangular faces where each face consists of three vertices.
%A simple form of the {\ct GEOM}\ namelist defining one triangular face is given by
%
%\begin{verbatim}
%&GEOM ID='triangle'
%      VERTS=0.0,0.0,0.0, 1.0,0.0,0.0, 1.0,0.0,1.0,
%      FACES=1,2,3
%/
%\end{verbatim}
%
%\noindent where {\ct ID}\ specifies the object name, in this case {\ct triangle},
%{\ct VERTS}\ specifies a list of one or more $(x,y,z)$ coordinates and {\ct FACES}\ specifies a list of vertices, 3 vertex indices for each
%face. Each index ranges from 1 to the number of vertices found on this {\ct GEOM}\ line.
%
%One may also define particular kinds of geometric objects,  blocks using {\ct XB}, spheres using {\ct SPHERE\_ORIGIN}\ and {\ct SPHERE\_RADIUS}
%and a 2D surface using {\ct ZVALS}.
%FDS generate vertices and faces to represent these objects, equivalent to what would have been defined using {\ct VERTS}\ and {\ct FACES}.
%
%\subsubsection{Blocks}
%A {\ct GEOM}\ namelist defining a block is given by
%
%\begin{verbatim}
%&GEOM ID='block'
%      XB=0.0,1.0,0.0,1.0,0.0,1.0 /
%\end{verbatim}
%
%\noindent where {\ct XB=xmin,xmax,ymin,ymax,zmin,zmax} defines the min and max bounds of the block.
%The {\ct XB}\ parameter is used in the same way as on an {\ct OBST}\ or {\ct VENT}\ line.
%A block may be discretized into many parts by
%specifying the {\ct IJK} parameter.  For example, {\ct IJK=8,6,4} would split the block
%into 8 parts along the $x$ dimension, 6 parts along the $y$ dimension and 4 parts along
%the $z$ dimension. By default, blocks are discretized
%so that the block faces are consistent in size with the grid resolution.
%
%\subsubsection{Spheres}
%A {\ct GEOM}\ namelist defining a sphere centered at $(0,0,0)$ with radius $1$ is given by
%
%\begin{verbatim}
%&GEOM ID='sphere'
%      SPHERE_ORIGIN=0.0,0.0,0.0,SPHERE_RADIUS=1.0 /
%\end{verbatim}
%
%\noindent Spheres are discretized by default so that each face is consistent in size with the grid resolution.
%One may specify a {\ct N\_LEVEL}\ parameter which specifies the number of times the sphere is split.
%
%\subsubsection{2D Surfaces}
%A {\ct GEOM}\ namelist defining a 2D surface is given by
%
%\begin{verbatim}
%&GEOM ID='terrain'
%      IJK=ivals,jvals,XB=xmin,xmax,ymin,ymax,ZVALS=..... /
%\end{verbatim}
%
%\noindent where {\ct XB}\ defines a rectangular region bounded by xmin, xmax, ymin and ymax where elevation data is defined with {\ct ZVALS}.
%{\ct IJK}\ specifies how many vertices occur in this region along each dimension in this region.
%In this example {\ct ivals}\ values occur along the x dimension and {\ct jval}\ values occur along the y dimension.
%The {\ct ZVALS} keyword is used to specify elevations at each $(x,y)$ location.
%The elevation data specified after the {\ct ZVALS}\ keyword is arranged in row major order.
%The first row contains {\ct ivals} elevation values occurring at the the ymax position from xmin to xmax.
%Similarly, the last row contains
%{\ct ivals}\ elevation values for the ymin position again from xmin to xmax.
%There are then {\ct jvals}\ rows and {\ct ivals}\ columns of elevation data.
%As with the blocks and spheres, FDS uses the information
%provided by these keywords to construct
%vertices and triangular faces.

%\subsection{Defining Solids}
%Unstructured geometric solids are also defined using the {\ct GEOM}\ namelist.
%These solids consist of a collection of tetrahedrons where each tetrahedron consists of four vertices.
%A simple form of the {\ct GEOM}\ namelist defining one tetrahedron is given by
%
%\begin{verbatim}
%&GEOM ID='tetrahedron'
%      VERTS=0.0,0.0,0.0, 1.0,0.0,0.0, 1.0,1.0,0.0,  1.0,1.0,1.0,
%      VOLUS=1,2,3,4
%/
%\end{verbatim}
%
%\noindent where {\ct ID}\ specifies the object name, in this case {\ct tetrahedron},
%{\ct VERTS}\ specifies a list of one or more $(x,y,z)$ coordinates and {\ct VOLUS}\ specifies a list of vertices, 4 vertex indices for each
%tetrahedron. Each index ranges from 1 to the number of vertices found on this {\ct GEOM}\ line.


\subsection{Transforming Objects}
An object may be transformed.
It  may be translated using {\ct XYZ}, scaled using {\ct SCALE} and rotated about {\ct XYZ0}\ using {\ct AZIM}\ and/or {\ct ELEV}.
For example, in the following {\ct GEOM}\ namelist,

{\small
\begin{verbatim}
&MOVE ID='ChairMove', AXIS=...,...,..., ROTATION_ANGLE=..., DX=, /
&GEOM ID='chair'
      SURF_ID='woodSID','apholsterySID',
      VERTS=X1,Y1,Z1,...,XM,YM,ZM,
      FACES=F1_1,F1_2,F1_3,F1_SID...,FN_1,FN_2,FN_3,F2_SID
      MOVE_ID='ChairMove'
/
\end{verbatim}
}

\noindent the chair object is translated and rotated using a reference to the $\&MOVE$ namelist, defined and explained in section XX.


\subsection{Reading Geometry node locations and connectivity data from binary}
\label{subsec:readbin}

In cases where geometries have very large numbers of nodes and surface triangles, it can be more efficient at setup time to read this data from an FDS produced binary file. This binary file can be obtained running the case in setup only mode ({\ct T\_END=0.0}).
Consider the chair geometry defined in the previous section. When a geometry is read by FDS, the {\ct XYZ} and face connectivity data {\ct FACES} is automatically written into a binary file with name {\ct CHID\_GEOM\_ID.bingeom} on the run directory (here, {\ct GEOM\_ID='chair'}).

Then, the next time the case is run, the geometry line can be modified to:
{\small
\begin{verbatim}
&GEOM ID='chair',
      SURF_ID='woodSID','apholsterySID',
      READ_BINARY=.TRUE.
      MOVE_ID='ChairMove'
/
\end{verbatim}
}

where {\ct READ\_BINARY=.TRUE.} instructs FDS to read nodes and faces from the binary geometry file.
Note we have deleted the {\ct VERTS, FACES} arrays from the input line.
This manner of loading large geometries avoids the testing and reading of all lines related to these arrays on the FDS input file.
It is particularly efficient in large parallel calculations, where parallel file systems are used. See section~\ref{inout:bingeom} for more information on geometry I/O binary format.


%\subsection{Grouping Objects}
%One may group objects to form a new object.  In the following example,
%{\ct chair}\ and {\ct couch}\ are objects that have been defined previously.  A new {\ct living room set}\ object
%is defined by referencing {\ct chair}\ and {\ct couch}\ in the {\ct GEOM\_IDS}\ keyword.  The chair and couch
%are place at a particular locations using the {\ct DXYZ}\ keyword.
%
%{\small
%\begin{verbatim}
%&GEOM ID='chair' ..... /
%
%&GEOM ID='sofa' ...... /
%
%&GEOM ID='living room set'
%
%      XYZ=[0,0,0],
%      AZIM=[0],
%      ELEV=[0],
%
%      GEOM_IDS(1)='chair',  DSCALE(1)=..., DAZIM(1)=..., DELEV(1)...
%                            DXYZ0(1:3,1)=..., DXYZ(1:3,1)=...,
%      GEOM_IDS(2)='couch',  DSCALE(2)=..., DAZIM(2)=..., DELEV(2)...
%                            DXYZ0(1:3,2)=..., DXYZ(1:3,2)=...,
%/
%\end{verbatim}
%}
%
%\noindent Referenced objects such as {\ct chair}\ and {\ct couch} in this case must be defined in
%previous {GEOM}\ lines.  Forward references are not permitted.
%One may scale, rotate and translate objects as they are placed in the group using
%{\ct DSCALE}\, for scaling, {\ct DAZIM} and {\ct DELEV} for rotation and {\ct DXYZ} for translation.
%
%The {\ct DXYZ0}\ and {\ct DXYZ}\ keywords use two coordinates (array indices) to specify data.
%The first coordinate specifies the spatial component (1 for x, 2 for y, 3 for z or 1:3 for all three).
%The second coordinate specifies the object number according to its position in the {\ct GEOM\_IDS}\ array.
%For example,
%DXYZ(1,3) would specify the translation along the x direction for the third object in the {\ct GEOM\_IDS}\ array.
%
%The {\ct COMPONENT\_ONLY}\ keyword if .TRUE., indicates that the object being defined is only a component found in other geometries.
%It will not appear by itself in the FDS model.

%\subsection{Dynamic Objects}
%An object's position and orientation may change with time.
%These changes are specified using various keywords that end with {\ct \_DOT}.
%See Table \ref{tbl:GEOM} for a complete list of these keywords.
%In the following example, {\ct AZIM\_DOT} is used to vary the azimuthal rotation angle by 1 deg/s.
%Similarly, the position may be changed by using {\ct XYZ\_DOT}.
%For example, {\ct XYZ\_DOT=0.0,0.0,1.0}\ would cause an object to move 1~m/s in a vertical direction.
%{\small
%\begin{verbatim}
%&GEOM VERTS=...,FACES=....,AZIM=0.0,AZIM_DOT=1.0 /
%\end{verbatim}
%}
%
%\vspace{\baselineskip}

\section{\texorpdfstring{{\ct GEOM}}{GEOM} (Unstructured Geometry Parameters)}

\begin{longtable}{@{\extracolsep{\fill}}|l|l|l|l|l|}
\caption[Unstructured geometry parameters ({\ct GEOM} namelist group)]{For more information see Section~\ref{info:GEOM}.}
\label{tbl:GEOM} \\
\hline
\multicolumn{5}{|c|}{{\ct GEOM} (Unstructured Geometry Parameters)} \\
\hline \hline
\endfirsthead
\caption[]{Continued} \\
\hline
\multicolumn{5}{|c|}{{\ct GEOM} (Unstructured Geometry Parameters)} \\
\hline \hline
\endhead
{\ct CYLINDER\_LENGTH}& Real                  & Section~\ref{info:GEOM}            &   m       &                      \\ \hline
{\ct CYLINDER\_RADIUS}& Real                  & Section~\ref{info:GEOM}            &   m       &                       \\ \hline
{\ct CYLINDER\_ORIGIN}& Real Triplet          & Section~\ref{info:GEOM}            &   m       &             \\ \hline
{\ct CYLINDER\_AXIS}& Real Triplet          & Section~\ref{info:GEOM}            &          &             \\ \hline
{\ct CYLINDER\_NSEG\_AXIS}& Integer          & Section~\ref{info:GEOM}            &          &             \\ \hline
{\ct CYLINDER\_NSEG\_THETA}& Integer          & Section~\ref{info:GEOM}            &          &             \\ \hline
{\ct EXTEND\_TERRAIN}& Logical          & Section~\ref{info:GEOM}            &          &    {\ct .FALSE.}         \\ \hline
{\ct FACES}        & array of Integer Triplets     & Section~\ref{info:GEOM}     &           &    0                     \\ \hline
{\ct ID}           & Character              & Section~\ref{info:GEOM}            &           &   {\ct 'geom'\_\#}           \\ \hline
{\ct IJK}          & Integer Triplet        & Section~\ref{info:GEOM}            &           &   0,0,0                  \\ \hline
{\ct MATL\_ID}     & Character              & Section~\ref{info:GEOM}            &           &  {\ct 'INERT'}           \\ \hline
{\ct MOVE\_ID}     & Character              & Section~\ref{info:GEOM}            &           &  {\ct 'null'}           \\ \hline
{\ct N\_LAT}       & Integer                & Section~\ref{info:GEOM}            &           &   0                      \\ \hline
{\ct N\_LEVELS}    & Integer                & Section~\ref{info:GEOM}            &           &   0                      \\ \hline
{\ct N\_LONG}      & Integer                & Section~\ref{info:GEOM}            &           &   0                      \\ \hline
{\ct READ\_BINARY}   & Logical         & Section~\ref{info:GEOM}            &              &   {\ct .FALSE.} \\ \hline
{\ct SPHERE\_ORIGIN}& Real Triplet          & Section~\ref{info:GEOM}            &   m       &             \\ \hline
{\ct SPHERE\_RADIUS}& Real                  & Section~\ref{info:GEOM}            &   m       &                    \\ \hline
{\ct SPHERE\_TYPE} & Integer                & Section~\ref{info:GEOM}            &           &                        \\ \hline
{\ct SURF\_ID}     & Character              & Section~\ref{info:GEOM}            &           &  {\ct 'INERT'}           \\ \hline
{\ct TEXTURE\_MAPPING}& Character           & Section~\ref{info:GEOM}            &           & {\ct 'RECTANGULAR'}      \\ \hline
{\ct TEXTURE\_ORIGIN} & Real Triplet        & Section~\ref{info:GEOM}            &   m       &   0.0,0.0,0.0            \\ \hline
{\ct TEXTURE\_SCALE}& Real                  & Section~\ref{info:GEOM}            &           &   1.0                    \\ \hline
{\ct VERTS}        & array of Real Triplets & Section~\ref{info:GEOM}            &   m       &   0.0                    \\ \hline
{\ct VOLUS}        & array of Integer Quadruplets     & Section~\ref{info:GEOM}  &           &    0                     \\ \hline
{\ct XB}           & Real sex-tuplet        & Section~\ref{info:GEOM}            &   m       &   0.0                    \\ \hline
{\ct ZVALS}        & Real                   & Section~\ref{info:GEOM}            &   m     &   0.0                    \\ \hline
{\ct ZVAL\_HORIZON}        & Real                   & Section~\ref{info:GEOM}    &      &                    \\ \hline


% Related to moving geometry:
%{\ct AZIM}         & Real                   & Section~\ref{info:GEOM}            &  deg      &    0.0                   \\ \hline
%{\ct COMPONENT\_ONLY} & LOGICAL             & Section~\ref{info:GEOM}            &           &  {\ct .FALSE.}           \\ \hline
%{\ct DAZIM}        & array of Reals         & Section~\ref{info:GEOM}            &  deg      &    0.0                   \\ \hline
%{\ct DELEV}        & array of Reals         & Section~\ref{info:GEOM}            &  deg      &    0.0                   \\ \hline
%{\ct DSCALE}       & array of Real Triplets & Section~\ref{info:GEOM}            &           &   1.0                    \\ \hline
%{\ct DXYZ0}        & array of Real Triplets & Section~\ref{info:GEOM}            &   m       &   0.0                    \\ \hline
%{\ct DXYZ}         & array of Real Triplets & Section~\ref{info:GEOM}            &   m       &   0.0                    \\ \hline
%{\ct ELEV}         & Real                   & Section~\ref{info:GEOM}            &  deg      &    0.0                   \\ \hline
%{\ct GAXIS}        & Real Triplet           & Section~\ref{info:GEOM}            &           &                          \\ \hline
%{\ct GEOM\_IDS}    & Real                   & Section~\ref{info:GEOM}            &           &                          \\ \hline
%{\ct GROTATE}      & Real Triplet           & Section~\ref{info:GEOM}            &  deg      &    0.0                   \\ \hline
%{\ct SCALE}        & Real Triplet           & Section~\ref{info:GEOM}            &           &   1.0,1.0,1.0            \\ \hline
%{\ct XYZ0}         & Real Triplet           & Section~\ref{info:GEOM}            &   m       &   0.0                    \\ \hline
%{\ct XYZ}          & Real Triplet           & Section~\ref{info:GEOM}            &   m       &   0.0                    \\ \hline

%{\ct AZIM\_DOT}    & Real                   & Section~\ref{info:GEOM}            &  deg/s    &    0.0                   \\ \hline
%{\ct ELEV\_DOT}    & Real                   & Section~\ref{info:GEOM}            &  deg/s    &    0.0                   \\ \hline
%{\ct GROTATE\_DOT}  & Real Triplet          & Section~\ref{info:GEOM}            & deg/s     &    0.0                   \\ \hline
%{\ct SCALE\_DOT}   & Real Triplet           & Section~\ref{info:GEOM}            &  1/s      &   1.0                    \\ \hline
%{\ct XYZ\_DOT}     & Real Triplet           & Section~\ref{info:GEOM}            &   m/s     &   0.0                    \\ \hline


\end{longtable}


\vspace{\baselineskip}

\section{File Formats}

\subsection{FDS GEOM I/O binary format}
\label{inout:bingeom}

The binary format for geometries read by FDS as described in section~\ref{subsec:readbin} is the following. Binary files read by the FDS released precompiled bundles must have \textit{little endian} endianness.
%
\begin{lstlisting}

WRITE(LU_BINGEOM) INTEGER_ONE
WRITE(LU_BINGEOM) N_VERTS,N_FACES,N_SURF_ID,N_VOLUS
WRITE(LU_BINGEOM) VERTS(1:3*N_VERTS)
WRITE(LU_BINGEOM) FACES(1:3*N_FACES)
WRITE(LU_BINGEOM) SURFS(1:N_FACES)
WRITE(LU_BINGEOM) VOLUS(1:4*N_VOLUS)

\end{lstlisting}

where:

\begin{itemize}
  \item {\ct INTEGER\_ONE} Integer equal to one, used to test that the endianness of the binary being read is \textit{little endian}.
  \item {\ct N\_VERTS,N\_FACES,N\_SURF\_ID,N\_VOLUS} number of vertices, triangles, surface types and volumes on the geometry.
  \item {\ct VERTS} one dimensional array with $x,y,z$ positions of the {\ct N\_VERTS} vertices.
  \item {\ct FACES} one dimensional array with the $n_1,n_2,n_3$ node connectivities of the {\ct N\_FACES} triangles.
  \item {\ct SURFS} one dimensional array with {\ct SURF\_ID}s of the {\ct N\_FACES} triangles.
  \item {\ct VOLUS} one dimensional array with $n_1,n_2,n_3,n_4$ node connectivities of the {\ct N\_VOLUS} tetrahedrons. This array is currently unused in FDS and left for 3D solid heat transfer in the future.
\end{itemize}


\subsection{SMV Unstructured Geometry}
\label{out:GEOMETRY}

Immersed geometric objects are stored for SMV reading using a file format described below.
These objects are defined in terms of vertices, triangles and tetrahedrons.
A vertex is represented as an $(x,y,z)$ coordinate (three floating point values).
A triangle is represented as three vertex (integer) indices. A tetrahedron is represented as four vertex (integer) indices.
An associated file format used to store data  on the geometric object is described in the next section.
These files are written out from {\ct dump.f90} using lines equivalent to the following:

\begin{lstlisting}
! header

WRITE(LU_GEOM) ONE
WRITE(LU_GEOM) VERSION
WRITE(LU_GEOM) N_FLOATS, N_INTS, FIRST_FRAME_STATIC
IF (N_FLOATS>0) WRITE(LU_GEOM) (FLOAT_HEADER(I),I=1,N_FLOATS)
IF (N_INTS>0) WRITE(LU_GEOM) (INT_HEADER(I),I=1,N_INTS)

! geometry frame
! STIME ignored if first frame is static ( FIRST_FRAME_STATIC set to 1)

WRITE(LU_GEOM) STIME
WRITE(LU_GEOM) N_VERT, N_FACE, N_VOL
IF (N_VERT>0) WRITE(LU_GEOM)(Xvert(I),Yvert(I),Zvert(I),I=1,N_VERT)
IF (N_FACE>0) THEN
   WRITE(LU_GEOM) (FACE1(I),FACE2(I),FACE3(I),I=1,N_FACE)
   WRITE(LU_GEOM) (SURF(I),I=1,N_FACE)
   WRITE(LU_GEOM) (Xtext(I),Ytext(I),I=1,3*N_FACE)
ENDIF
IF (N_VOL>0) THEN
   WRITE(LU_GEOM) (VOL1(I),VOL2(I),VOL3(I),VOL4(I),I=1,N_VOL)
   WRITE(LU_GEOM) (MATL(I),I=1,N_VOL)
ENDIF
              .
\end{lstlisting}

\begin{itemize}
\item {\ct ONE}\ has the value 1. Smokeview uses this number to determine whether the computer creating the geometry file and the computer viewing the geometry file use the same or different byte swap (endian) conventions for storing floating point numbers.
\item {\ct VERSION}\ currently has value 2 and indicates the version number of this file format.
\item {\ct N\_FLOATS, N\_INTS} The number of floating point and integer data items stored at the beginning of the file.
\item {\ct FLOAT\_HEADER, INT\_HEADER} Floating point and integer data stored at the beginning of the file.
\item {\ct STIME} is the FDS simulation time.
\item {\ct N\_VERT, N\_FACE, N\_VOL}  are the number of vertices, faces and volumes.
\item {\ct Xtext, Ytext}\ are the texture coordinates.
\item {\ct Xvert, Yvert, Zvert}\ are the vertex coordinates.
\item {\ct FACE1, FACE2, FACE3}\ are the vertex indices for each face (triangle).
    The indices range from 1 to the
    number of vertices.
\item {\ct VOL1, VOL2, VOL3, VOL4}\ are the vertex indices for each volume (tetrahedron).  The indices are numbered from 1 to the number of vertices.
\item {\ct SURF}\ are the SURF indices for each face.
\item {\ct MATL}\ are the MATL indices for each volume.
\end{itemize}


\subsection{Unstructured Data}

\label{out:GEOMETRYDATA}

(add to FDS user's guide 20.13 when ready)

This section describes the file format used for
storing data associated with unstructured geometric objects.  This format will be used to unify file formats for storing 3D smoke, boundary, particle and slice data.
Unstructured data files are written out from {\ct dump.f90} using lines equivalent to the following:

\begin{lstlisting}

! header

WRITE(LU_GEOM) ONE
WRITE(LU_GEOM) VERSION
WRITE(LU_GEOM) N_FLOATS
IF (N_FLOATS>0) WRITE(LU_GEOM) (FLOAT_HEADER(I),I=1,N_FLOATS)
WRITE(LU_GEOM) N_INTS
IF (N_INTS>0) WRITE(LU_GEOM) (INT_HEADER(I),I=1,N_INTS)

! data for each time step

WRITE(LU_GEOM) STIME
WRITE(LU_GEOM) N_VERT, N_FACE, N_VOL
WRITE(LU_GEOM) HAS_VERT_ID
IF (N_VERT>0) THEN
   WRITE(LU_GEOM)(ValVert(I),I=1,N_VERT)
   IF (HAS_VERT_ID.EQ.1) WRITE(LU_GEOM)(Vert_ID(I),I=1,N_VERT)
ENDIF
IF (N_FACE>0) WRITE(LU_GEOM)(ValFace(I),I=1,N_FACE)
IF (N_VOLS>0) WRITE(LU_GEOM)(ValVol(I),I=1,N_VOLS)

\end{lstlisting}

\begin{itemize}
\item {\ct ONE}\ has the value 1. Smokeview uses this number to determine whether the computers creating and viewing this data file use the same or different conventions for storing floating point numbers.
\item {\ct VERSION}\ version number of the format used to store data in this file.
\item {\ct N\_FLOATS, N\_INTS} The number of floating point and integer data items stored at the beginning of the file.
\item {\ct FLOAT\_HEADER, INT\_HEADER} Floating point and integer data stored at the beginning of the file.
\item {\ct STIME} is the simulation time.
\item {\ct N\_VERT, N\_FACE, N\_VOLS}\ are the number of data values associated with vertices, faces and volumes.
\item {\ct ValVert}\ is data associated with the vertices in the unstructured geometric object.
\item {\ct ValFace}\ is data associated with the faces in the unstructured geometric object.
\item {\ct ValVol}\ is data associated with the volumes in the unstructured geometric object.
\end{itemize}

\subsection{Isosurface Data}
\label{out:ISOSURFACE}

FDS generated isosurfaces are stored using a
file format described in this section.
Iso-surface files are used to store one or more surfaces where the
surface has the value specified in the {\ct QUANTITY}\ keyword.
FDS outputs iso-surface data at fixed time intervals.
These surfaces are defined in terms of vertices and triangles.
A vertex is represented as an $(x,y,z)$ coordinate (three floating point values).
A triangle is represented as three vertex (integer) indices.
These files are written out from {\ct dump.f90} using lines equivalent to the following:
\begin{lstlisting}
! header

WRITE(LU_GEOM) ONE
WRITE(LU_GEOM) VERSION
WRITE(LU_GEOM) N_FLOATS
IF (N_FLOATS>0) WRITE(LU_GEOM) (FLOAT_HEADER(I),I=1,N_FLOATS)
WRITE(LU_GEOM) N_INTS
IF (N_INTS>0) WRITE(LU_GEOM) (INT_HEADER(I),I=1,N_INTS)

! static geometry - geometry specified once and appearing at all time steps

WRITE(LU_GEOM) N_VERT_S, N_FACE_S
IF (N_VERT_S>0)  WRITE(LU_GEOM) (Xvert_S(I),Yvert_S(I),Zvert_S(I),I=1,N_VERT_S)
IF (N_FACE_S>0)  THEN
   WRITE(LU_GEOM) (FACE1_S(I),FACE2_S(I),FACE3_S(I),I=1,N_FACE_S)
   WRITE(LU_GEOM) (SURF_S(I),I=1,N_FACE_S
ENDIF

! dynamic geometry - geometry specified and appearing for each time step

WRITE(LU_GEOM) STIME, unused
WRITE(LU_GEOM) N_VERT_D, N_FACE_D
IF (N_VERT_D>0) WRITE(LU_GEOM)(Xvert_D(I),Yvert_D(I),Zvert_D(I),I=1,N_VERT_D)
IF (N_FACE_D>0) THEN
   WRITE(LU_GEOM) (FACE1_D(I),FACE2_D(I),FACE3_D(I),I=1,N_FACE_D)
   WRITE(LU_GEOM) (SURF_D(I),I=1,N_FACE_D)
ENDIF
\end{lstlisting}

\begin{itemize}
\item {\ct ONE}\ has the value 1. Smokeview uses this number to determine whether the computer creating the geometry file and the computer viewing the geometry file use the same or different byte swap (endian) conventions for storing floating point numbers.
\item {\ct VERSION}\ has value 0 (for isosurface files but has value 1 for geometry files) and indicates the version number of this file format.
\item {\ct N\_FLOATS, N\_INTS} The number of floating point and integer data items stored at the beginning of the file.
\item {\ct FLOAT\_HEADER, INT\_HEADER} Floating point and integer data stored at the beginning of the file.
\item {\ct STIME} is the FDS simulation time.
\item {\ct N\_VERT\_S, N\_FACE\_S, N\_VERT\_D, N\_FACE\_D} are the number of static and dynamic vertices and faces.
\item {\ct Xvert\_S, Yvert\_S, Zvert\_S, Xvert\_D, Yvert\_D, Zvert\_D}\ are the static and dynamic vertex coordinates.
\item {\ct FACE1\_S, FACE2\_S, FACE3\_S, FACE1\_S, FACE2\_S, FACE3\_S}\ are the static and dynamic vertex indices for each face (triangle).
    The indices range from 1 to the number of vertices.
\item {\ct SURF\_S, SURF\_D}\ are the static and dynamic SURF indices for each face (color indices for isosurface files.
\end{itemize}









\section{Computational Geometry Technical Notes:}

The notes written in this section are intended to provide details for various computational algorithms as implemented in FDS and Smokeview.


%
%\subsection{Geometry Properties}
%
%
%
%\subsubsection{Triangle,Polygon}
%Consider the triangle with vertices $v_1$,$v_2$ and $v_3$ as illustrated in Figure \ref{figure:triangle_setup}.
%
%\noindent Triangle area: $||(v_3-v_1)\times (v_2-v_1)||/2$ . \\
%Triangle Median: $(v_1+v_2+v_3)$/3 .
%
%\begin{figure}
%\begin{center}
%\includegraphics[width=4.0in]{../../../fig/smv/figures/triangle_setup}
%\end{center}
%\caption{Schematic of vertex notation for a triangle}
%\label{figure:triangle_setup}
%\end{figure}
%
%\subsubsection{Tetrahedron,Polyhedron}
%Consider the tetrahedron T with vertices $v_1$,$v_2$,$v_3$ and $v_4$ as illustrated in Figure \ref{figure:tetrahedron_setup}.
%
%\noindent Tetrahedron volume: $||((v_3-v_1)\times (v_2-v_1))\cdot (v_4-v_1)||/6$ .\\
%Tetrahedron median: $(v_1+v_2+v_3+v_4)$/4.
%
%\begin{figure}
%\begin{center}
%\includegraphics[width=4.0in]{../../../fig/smv/figures//tetrahedron_setup}
%\end{center}
%\caption{Schematic of vertex and edge notation for a tetrahedron}
%\label{figure:tetrahedron_setup}
%\end{figure}
%
%
%\newcommand{\tetra}{\mbox{tetra\_bounds}}
%\newcommand{\bbox}{\mbox{box\_bounds}}
%\newcommand{\txtmin}{\mbox{min}}
%\newcommand{\txtmax}{\mbox{max}}
%\subsection{Tetrahedron Box Volume Intersection}
%\subsubsection{Overview}
%The region formed by intersecting a box and a tetrahedron may have up to 10 faces.  Up to four faces occur on the tetrahedron and up to six faces occur on the box.  Each box (or tetrahedron) face is a part of the intersection region if it intersects with the tetrahedron (or box).  The problem then is to determine the volume and surface area of this intersection region.  A general algorithm for computing the intersection volume is given below with more details given in subsequent sections.
%\begin{enumerate}
%\item Determine the tetrahedron and box bounding boxes.  If the bounding boxes do not overlap then the intersection volume is zero.  Note, if the bounding boxes do overlap the intersection volume may still be zero.
%
%\item Determine vertices in intersection region occurring on each box face $b_i$:
%\label{boxstep}
%\begin{enumerate}
%\item For each tetrahedron face $t_j$:
%\label{tetraverts}
%\begin{enumerate}
%\item Intersect $b_i$ with each $t_j$ edge.
%\item Add any vertices of $b_i$ that are on or inside the tetrahedron.
%\end{enumerate}
%
%\item If two or fewer vertices are found in step \ref{tetraverts} then ignore this face since
%the volume of this portion of the intersection region is zero.  This may occur if the box and tetrahedron  intersect along an edge.
%\item If three or more vertices are found then order them clockwise.
%\end{enumerate}
%
%
%\item Determine vertices in intersection region occurring on each tetrahedron face $t_i$:
%\label{tetrastep}
%\begin{enumerate}
%\item For each box face $b_j$ not in same plane as $t_i$:
%\label{boxverts}
%
%\begin{enumerate}
%\item Intersect $t_i$ with each $b_j$ edge.
%\item Add any $t_i$ vertices on or inside the box.
%\end{enumerate}
%Note, vertices for this face were already found in step \ref{tetraverts} if $b_j$ and $t_i$ are in the same plane.
%\item If two or fewer vertices are found in step \ref{boxverts} then ignore this face since the volume of this portion of the intersection region is zero. As before, this may occur if the box and tetrahedron  intersect along an edge.
%\item If three or more vertices are found then order them clockwise.
%\end{enumerate}
%
%\item Form an interior vertex, $v_\mathrm{int}$, by averaging all vertices found in steps \ref{boxstep} and \ref{tetrastep}.
%
%\item For each intersection region face $f_i$, compute the volume formed by appending $f_i$ with $v_\mathrm{int}$.
%\label{volfacestep}
%
%\item Sum the volumes computed for each face $f_i$ in step \ref{volfacestep}.
%\end{enumerate}
%
%\subsubsection{Determine Bounding Boxes}
%\begin{enumerate}
%\item Bounding box for a tetrahedron (see Figure \ref{figure:tetrahedron_setup} for box vertex definitions):
%\begin{eqnarray*}
%\tetra(0)&=x_{\txtmin}=\min(v0_x,v1_x,v2_x,v3_x)\\
%\tetra(1)&=x_{\txtmax}=\max(v0_x,v1_x,v2_x,v3_x)\\
%\tetra(2)&=y_{\txtmin}=\min(v0_y,v1_y,v2_y,v3_y)\\
%\tetra(3)&=y_{\txtmax}=\max(v0_y,v1_y,v2_y,v3_y)\\
%\tetra(4)&=z_{\txtmin}=\min(v0_z,v1_z,v2_z,v3_z)\\
%\tetra(5)&=z_{\txtmax}=\max(v0_z,v1_z,v2_z,v3_z)
%\end{eqnarray*}
%
%\item Bounding box for for an axis aligned box (see Figure \ref{figure:box_setup} for tetrahedron vertex definitions):
%\begin{eqnarray*}
%\bbox(0)&=x_{\txtmin}=v0_x\\
%\bbox(1)&=x_{\txtmax}=v1_x\\
%\bbox(2)&=y_{\txtmin}=v0_y\\
%\bbox(3)&=y_{\txtmax}=v2_y\\
%\bbox(4)&=z_{\txtmin}=v0_z\\
%\bbox(5)&=z_{\txtmax}=v4_z
%\end{eqnarray*}
%\end{enumerate}
%
%\noindent The intersection volume is zero if bounding boxes for the tetrahedron and box do not overlap.
%These bounding boxes do not over lap if any of the following are true:
%\begin{eqnarray*}
%\mbox{box\_bounds}(0)>\mbox{tetra\_bounds}(1) \\
%\mbox{box\_bounds}(1)<\mbox{tetra\_bounds}(0) \\
%\mbox{box\_bounds}(2)>\mbox{tetra\_bounds}(3) \\
%\mbox{box\_bounds}(3)<\mbox{tetra\_bounds}(2) \\
%\mbox{box\_bounds}(4)>\mbox{tetra\_bounds}(5) \\
%\mbox{box\_bounds}(5)<\mbox{tetra\_bounds}(4)
%\end{eqnarray*}
%
%\begin{figure}
%\begin{center}
%\includegraphics[width=4.0in]{../../../fig/smv/figures//box_setup}
%\end{center}
%\caption{Schematic of vertex and edge notation for a box}
%\label{figure:box_setup}
%\end{figure}
%
%\subsubsection{Distance from a Vertex to a Plane}
%\subsubsection{Intersection of a Plane with a Line Segment}
%\subsubsection{Determine if a Vertex is Inside a Box or Tetrahedron}
%\subsubsection{Intersection of a Rectangle and a Triangle Occurring in the Same Plane}
%\subsubsection{Orienting Vertices Occurring in the Same Plane}
%\subsubsection{Computing Volume of a Tetrahedron}
%
%\subsection{Spheres}
%\subsubsection{Latitude/Longitude Discretization}
%\subsubsection{Recursive Discretization}
%Vertex locations for initial discretization
%\begin{eqnarray*}
%v_1&=&(0,0,1)\\
%v_i&=&\left(\cos\left(2\pi\frac{i-1}{5}\right)\sqrt{5/6},\sin\left(2\pi\frac{i-1}{5}\right)\sqrt{5/6},\sqrt{1/6}\right)i=2,\cdots, 6\\
%v_i&=&\left(\cos\left(2\pi\frac{i-7}{5}\right)\sqrt{5/6},\sin\left(2\pi\frac{i-7}{5}\right)\sqrt{5/6},-\sqrt{1/6}\right)i=7,\cdots, 11\\
%v_{12}&=&(0,0,-1)
%\end{eqnarray*}
%
%Euler formula for relating the number vertices (V), edges (E) and faces (F) of a closed geometric surface used to discretize the sphere.
%\begin{eqnarray*}
%V+F-E&=&2
%\end{eqnarray*}
%Number of vertices, face and edges for the initial discretization.
%\begin{eqnarray*}
%V_0&=&12\\
%F_0&=&20\\
%E_0&=&30\\
%\end{eqnarray*}
%Recursion formulas for relating the number of vertices, faces and edges at the $n+1$'st refinement step given the corresponding number at the $n$'th refinement step.
%\begin{eqnarray}
%\label{eq:recursion1}
%V_{n+1}&=&V_n+E_n\\
%\label{eq:recursion2}
%F_{n+1}&=&4 F_n\\
%\label{eq:recursion3}
%E_{n+1}&=&3 F_n+2 E_n
%\end{eqnarray}
%
%For $n\ge 0$ show that $V_{n+1}+F_{n+1}-E_{n+1}=2$.
%Substituting $n=0$ into Eqs. (\ref{eq:recursion1}) through (\ref{eq:recursion3})  results in
%\begin{eqnarray*}
%V_1+F_1-E_1&=&V_0+E_0+4 F_0 -(3 F_0+2 E_0)\\
%&=&12+30+80-(60+60)\\
%&=&2
%\end{eqnarray*}
%
%Assuming the recursion formulas in Eqs. (\ref{eq:recursion1}) through (\ref{eq:recursion3}) and that $V_n+F_n-E_n=2$ (the induction step) results in
%\begin{eqnarray*}
%V_{n+1}+F_{n+1}-E_{n+1}&=&V_n+E_n+4 F_n-(3 F_n+2 E_n)\\
%&=&V_n+F_n-E_n\\
%&=&2
%\end{eqnarray*}
%
%Table \ref{table:refinement} gives the number of vertices, faces and edges for the 0'th through 5'th refinements.
%This table will be used to verify that the sphere discretization routine, {\ct INIT\_SPHERE}, found in FDS generates the correct number of faces and vertices..
%\vskip 0.1in
%\par
%\begin{table}
%\caption{Number of vertices, faces and edges for the 0'th through 5'th discretization refinements}
%\begin{tabular}{|c|c|c|c|}
%  \hline
%n& V&F&E\\
%\hline
%0&	12&	20&	30\\
%1&	42&	80&	120\\
%2&	162&	320&	480\\
%3&	642&	1280&	1920\\
%4&	2562&	5120&	7680\\
%5&	10242&	20480&	30720\\
%  \hline
%\end{tabular}
%\label{table:refinement}
%\end{table}
%
%Table \ref{table:facelist} gives the list of vertices for each face in
%in the initial refinement where the vertices are numbered as given  in
%Figure \ref{figure:facelist}.
%\begin{table}
%\caption{Vertex list for each face}
%\begin{center}
%\begin{tabular}{ccccc}
%  \hline
%  (1,2,3)&   (1, 3,4)&  (1, 4, 5)&   (1, 5, 6)&   (1, 6,2)\\
%  (2,11,7)&  (2,7,3)&   (3, 7,8)&  (3, 8, 4)&   (4, 8, 9)\\
%  (4, 9,5)&  (5,9,10)&   (5,10,6)& (6,10,11)&   (6,11, 2)   \\
%  (11,12,7)& (7,12,8)&   (8,12,9)&   (9,12,10)&   (10,12,11)   \\
%  \hline
%\end{tabular}
%\end{center}
%\label{table:facelist}
%\end{table}
%
%\begin{figure}
%\begin{center}
%\includegraphics[width=4.0in]{../../../fig/smv/figures/icosahedron_setup}
%\end{center}
%\caption{Schematic of vertex locations for defining face lists}
%\label{figure:facelist}
%\end{figure}
%
%\subsection{Boxes}
%\subsection{2D Surfaces}


\clearpage

\section{Cut Cell Immersed Boundary Method (CC\_IBM) Technical Notes}


\subsection{Introduction}


\subsection{Computational Geometry engine for cut-cell definition in FDS}


\subsubsection{Cut-cell definition scheme implementation notes}


\subsection{Scalar Transport discretization for Complex geometry}


\subsubsection{Multispecies mass balance equations}

Consider a set of gaseous, reacting  chemical species $\alpha=1,\dots,N$ flowing on a given spatial domain $\Omega \in \mathbb{R}^n, \; n=2,3$, with boundary $\partial \Omega$, paramaterized by an Eulerian reference frame $N$. These species are transported on a given point $\mathbf{x}$ in space with velocity $\mathbf{v}_\alpha(\mathbf{x},t)$ respect to $N$, and a mass weighted average velocity $\mathbf{v}(\mathbf{x},t)$
%
\begin{equation}
  \mathbf{v} = \frac{ \sum\limits_{\alpha=1}^{N} {\rho_\alpha \mathbf{v}_\alpha}}{\rho} \; , \; \rho =  \sum\limits_{\alpha=1}^{N} {\rho_\alpha} \label{eq:veldens}
\end{equation}
%
where space and time dependencies are not shown for simplicity, $\rho_\alpha(\mathbf{x},t) = \rho(\mathbf{x},t) Y_\alpha (\mathbf{x},t)$, $\rho$ is the mixture density and $Y_\alpha = \rho_\alpha / \rho$ is species $\alpha$ mass fraction. Definition: $\rho_\alpha(\mathbf{x},t)$ is the amount of mass of species $\alpha$ on a given differential volume $d\Omega$ centered in $\mathbf{x}$, while  $\rho(\mathbf{x},t)$ is the mass of \textit{all} species in $d\Omega$.

As defined, each species $\alpha$ has a velocity $\mathbf{v}_\alpha(\mathbf{x},t)$ which is different from the average mixture velocity $\mathbf{v}(\mathbf{x},t)$. This difference quantity is called the \textit{diffusion velocity} of species $\alpha$, $\mathbf{V}_\alpha$
%
\begin{equation}
   \mathbf{V}_\alpha(\mathbf{x},t) = \mathbf{v}_\alpha(\mathbf{x},t) - \mathbf{v}(\mathbf{x},t) \label{eq:vdiff}
\end{equation}
%

The transport-reaction mass balance equation for each individual species $\alpha$ is given by
%
\begin{equation}
   \frac{\partial \rho_\alpha}{ \partial t} + \nabla \cdot (\rho_\alpha  \mathbf{v}_\alpha) = \dot{m}_\alpha''' \; , \; \alpha=1,\dots,N \label{eq:bal}
\end{equation}
%
in $\Omega$, where $\dot{m}_\alpha'''(\mathbf{x},t)$ is the chemical reaction volume source (area source for $n=2$, sink if negative) of species $\alpha$. The problem is completely defined stating initial and boundary conditions for $\rho_\alpha(\mathbf{x},t)$ as well as time dependent fields $\mathbf{v}_\alpha(\mathbf{x},t)$, $\dot{m}_\alpha'''(\mathbf{x},t)$.

Noting that $ \rho_\alpha = \rho Y_\alpha$, and using equation~\eqref{eq:vdiff} in the previous
%
\begin{equation}
   \frac{\partial \rho Y_\alpha}{ \partial t} + \nabla \cdot \left( \rho Y_\alpha  (\mathbf{v}+\mathbf{V}_\alpha) \right) = \dot{m}_\alpha''' \; , \; \alpha=1,\dots,N \label{eq:bal2}
\end{equation}
%

We see from these equations that the evolution of these species can be considered as a function of the mixture $\rho$ and average $\mathbf{v}$, and written as such, depending on the other species through equations~\eqref{eq:veldens}, and through the mass conservation statement on reactive sources
%
\begin{equation}
  \sum\limits_{\alpha=1}^{N} \dot{m}_\alpha'''(\mathbf{x},t) = 0
\end{equation}
%

In equation~\eqref{eq:bal2}, the following convective and diffusive fluxes per unit time and area are found:
%
\begin{eqnarray}
  \mathbf{J_{c \alpha}} &=& \rho Y_\alpha  \mathbf{v} \label{eq:jc} \\
  \mathbf{J_{d \alpha}} &=& \rho Y_\alpha  \mathbf{V}_\alpha \label{eq:jd}
\end{eqnarray}
%
The definition of diffusion velocities $\mathbf{V}_\alpha$ is a central component of multispecies mass transport. For atmospheric combustion we consider the binary diffusion of all species respect to Nitrogen, the most abundant (\textit{background}) species. This simplification allows us to uncouple the diffusive fluxes, using Ficks Law of diffusion
%
\begin{equation}
   \mathbf{J_{d \alpha}} = \rho Y_\alpha  \mathbf{V}_\alpha = - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \label{eq:fick}
\end{equation}
%
where $D_\alpha$ is the binary diffusivity coefficient respect to the background species.

Using this last expression in equation~\eqref{eq:bal2} we arrive to the basic form of the species balance equations based on mass fractions, and used in FDS
%
\begin{equation}
   \frac{\partial \rho Y_\alpha}{ \partial t} + \nabla \cdot ( \rho Y_\alpha  \mathbf{v} ) = \nabla \cdot ( \rho D_\alpha \boldsymbol{\nabla} Y_\alpha ) + \dot{m}_\alpha''' \; , \; \alpha=1,\dots,N \label{eq:bal3}
\end{equation}
%

A spatially discretized version of the previous equations in a domain divided in $n_{tot}$ computational cells is given by:

%
\begin{equation}
\left[ \mathbf{M} \right] \frac{\partial}{\partial t} \left\{ \mathbf{\rho Y_\alpha} \right\} + \left[ \mathbf{A}_{adv}  \right] \left\{ \mathbf{\rho Y_\alpha} \right\} = \left\{ \mathbf{F_{diff \alpha}} \right\} + \left\{  \mathbf{\dot{m}_\alpha'''} \right\} +  \left\{ \mathbf{F_{BC \alpha}} \right\} \; , \; \alpha=1,\dots,N \label{eq:discbal3}
\end{equation}
%
where the vector $\left\{ \mathbf{\rho Y_\alpha} \right\}$ of size $n_{tot}$ is the vector of unknowns (assuming one unknown per computational cell), matrices $\left[ \mathbf{M} \right]$ and $ \left[ \mathbf{A_{adv}}  \right]$ are the systems mass and advection matrices, and $\left\{ \mathbf{F_{diff \alpha}} \right\}$, $\left\{  \mathbf{\dot{m}_\alpha'''} \right\}$, $\left\{ \mathbf{F_{BC \alpha}} \right\}$ are vectors due to diffusion, reaction and boundary conditions. Note that $\left\{ \mathbf{F_{diff \alpha}} \right\}$ is a function of the equations unknowns $\rho Y_\alpha$. In the next section we recast this term in a format amenable to be defined discretely as matrix-vector products.


\subsubsection*{Diffusive flux decomposition for Implicit time integration}

Consider the exact decomposition of the diffusive flux in equations~\eqref{eq:bal3}
%
\begin{equation}
   \mathbf{J_{d \alpha}} = - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha =  - \left( D_\alpha \boldsymbol{\nabla} ( \rho Y_\alpha)
   - \frac{D_\alpha}{\rho} \boldsymbol{\nabla} \rho \; ( \rho Y_\alpha) \right) \label{eq:expdfl2}
\end{equation}
%
which, substituted in the balance equation~\eqref{eq:bal3} leads to:
%
\begin{equation}
 \frac{\partial \rho Y_\alpha}{ \partial t} + \nabla \cdot \left(  \mathbf{v}' (\rho Y_\alpha)   \right) = - \nabla \cdot \left(  \mathbf{J_{d \alpha}'} \right) + \dot{m}_\alpha''' \; , \; \alpha=1,\dots,N \label{eq:bal4}
\end{equation}
%
where
%
\begin{eqnarray}
  \mathbf{v}' &=& \mathbf{v} + \frac{D_\alpha}{\rho} \boldsymbol{\nabla} \rho  \label{eq:vprime} \\
   \mathbf{J_{d \alpha}'} &=& -D_\alpha \boldsymbol{\nabla}  (\rho Y_\alpha )  \label{eq:jdaprime}
\end{eqnarray}
%
are the sum $\mathbf{v}'$ of the fluid mixture velocity $\mathbf{v}$ and a diffusion driven velocity term for species $\alpha$ in presence of density gradients $D_\alpha \boldsymbol{\nabla} \rho / \rho $, and a constant density diffusive flux term $\mathbf{J_{d \alpha}'}$.

A spatially discretized version of equations~\eqref{eq:bal4} in $n_{tot}$ computational cells is given by:
%
\begin{equation}
\left[ \mathbf{M} \right] \frac{\partial}{\partial t} \left\{ \mathbf{\rho Y_\alpha} \right\} + \left( \left[ \mathbf{A}_{adv \alpha}'  \right] +\left[ \mathbf{A}_{diff \alpha}' \right] \right) \left\{ \mathbf{\rho Y_\alpha} \right\}   =  \left\{  \mathbf{\dot{m}_\alpha'''} \right\} +  \left\{ \mathbf{F_{BC \alpha}} \right\} \; , \; \alpha=1,\dots,N \label{eq:discbal4}
\end{equation}
%
where in this case $ \left[ \mathbf{A}_{adv \alpha}'  \right]$ is the advection matrix due to the velocity field $\mathbf{v}'$, and  $\left[ \mathbf{A}_{diff \alpha}' \right]$ is the diffusion matrix due to $ \mathbf{J_{d \alpha}'} $. Note also that boundary conditions might not only produce right hand side entries, but also affect the advection and diffusion discretization matrices.



%The spatially discretized version of system~\eqref{}



\subsubsection{Finite Volume method on staggered grid in 2D}

The finite volume discretization (FV) starts by considering the integral form of equation~\eqref{eq:bal3} over a cell control volume $\Omega_{ii}$. For
cell $ii$ individualized by the index pair $(i,j)$ we have
%
\begin{equation}
 \int_{\Omega_{ii}} {\frac{\partial \rho Y_\alpha}{\partial t}} d \Omega + \int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot  \left(  \rho Y_\alpha \mathbf{v} \right)
      } d \Omega  = -\int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot \left(  \mathbf{J_{d \alpha}}  \right)  } d \Omega + \int_{\Omega_{ii}} { \dot{m}_\alpha''' } d \Omega \label{eq:intconvdiff}
\end{equation}
%
Assuming a time independent control volume, the time derivative and source terms are approximated by
%
\begin{eqnarray}
  \int_{\Omega_{ii}} {\frac{\partial \rho Y_\alpha}{\partial t}} d \Omega & = & \frac{\partial}{\partial t} \int_{\Omega_{ii}} {\rho Y_\alpha} d \Omega
  = \frac{\partial \widetilde{\rho \: Y_\alpha }_{i,j}}{\partial t} V_{i,j} \\
  \int_{\Omega_{ii}} { \dot{m}_\alpha''' } d \Omega & = & \widetilde{ \dot{m}_\alpha''' }_{i,j} V_{i,j} \label{eq:intcons}
\end{eqnarray}
%
where $V_{i,j}$ is the volume of cell $(i,j)$. These cell averaged quantities (denoted by tildes) match the cell centroid values of the corresponding scalar fields up to second order spatial accuracy. A consequence of this known fact is that using a second order finite difference method (FD), or a \textit{difference} finite volume approach with same definition
 for diffusion flux spatial derivatives and interpolation on the advective term, will lead to identical discretization matrices (they will vary up to a volume factor) on uniform cartesian grids. The difference between these two numerical discretization methods arises from the source $\dot{m}_\alpha''' $, which for a given cell is the cell centroid value $\dot{m}'''_{\alpha (i,j)}$ in FD, and it is the cell average $ \widetilde{ \dot{m}_\alpha''' }_{i,j}$ in FV. To discretize both diffusive and advective terms in equation~\eqref{eq:intconvdiff}, we make use of the divergence theorem. In the following we drop the tildes to simplify the notation, keeping in mind that using FV, quantities will always be cell or face averaged where it corresponds.

\subsubsection*{FV spatial discretization of Diffusive term} \label{Sec:FVdiff}

\subsubsection*{Ficks Law of diffusion}
Consider Fick's law governing the diffusive flux of quantity $\rho Y_\alpha$, $ \mathbf{J_{d \alpha}} = - \rho \: D_\alpha \boldsymbol{\nabla} Y_\alpha$. This expression is used in explicit integration of species transport. For cell $ii$, the first right hand side term of equation~\eqref{eq:intconvdiff} is
%
\begin{equation}
  \int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot \left(  \mathbf{J_{d \alpha}} \right)  } d \Omega =
  \int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)  } d \Omega =
  \int_{\partial \Omega_{ii}} { \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right) \cdot \hat{\mathbf{n}}_{ii} } \: d \partial \Omega
\end{equation}
%
where $\partial \Omega_{ii}$ is the boundary surface of cell $ii$, and $\hat{\mathbf{n}}_{ii}$ is the normal unit vector pointing outside of the cell control volume. Assuming a number of \textit{planar} faces $k=1,...,nf_c$ (with constant  $\hat{\mathbf{n}}_{ii,k}$) describe the surface of the cell, we have
%
\begin{equation}
    \int_{\partial \Omega_{ii}} { \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right) \cdot \hat{\mathbf{n}}_{ii} } \: d \partial \Omega = \sum^{nf_c}_{k=1}
    \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k \label{eq:discfvdiff}
\end{equation}
%
where the last term in parenthesis is the average value of $ \mathbf{J_{d \alpha}} $ over face $k$, and $A_k$ is the corresponding area.
%
\begin{figure}[h]
      %\centering
      \includegraphics[trim = 60mm 40mm 60mm 30mm, clip,
       width=0.49\linewidth]{../../../fig/fds/Geometry_Figures/CellsFVDiscretization.png}
      \includegraphics[trim = 60mm 40mm 60mm 30mm, clip,
       width=0.49\linewidth]{../../../fig/fds/Geometry_Figures/CellsFVHiXBc.png}
      \put(-350,-10){(a)}
      \put(-120,-10){(b)}
      \caption{Finite volume discretization in 2D: (a) Control volume of staggered cell $i,j$ individualized also by single index $ii$.  (b) Low $x$ boundary cell. In face $(\mathbf{ILO_f},j)$, corresponding to $\partial \Omega$, boundary conditions are derived.}
	\label{Fig:FVdisc}
\end{figure}
%
We use centered differences to estimate the $\boldsymbol{\nabla} Y_\alpha$ at each cell face. See figure~\ref{Fig:FVdisc}a. Being $[\rho D_\alpha]_{i,j}$ a cell centered quantity on the staggered cell, we have for the east face ($ii,k_e \rightarrow i+1/2,j$) of a generic 2D cell:
%
\begin{eqnarray}
   \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_{i+1/2,j} \cdot \hat{\mathbf{n}}_{ii,ke} \: A_{i+1/2,j} & = &
   -[\rho D_\alpha]_{i+1/2,j} \frac{Y_{\alpha i+1,j}-Y_{ \alpha i,j}}{\Delta x_{f(i)}} \Delta y_{c(j)} \label{eq:diffii}
\end{eqnarray}
%
For the other adjacent cell $jj$, with index pair $(i+1,j)$ the normal unit vector is $\hat{\mathbf{n}}_{jj,kw}=-\hat{\mathbf{n}}_{ii,ke}$. Therefore the dot product becomes:
%
\begin{eqnarray}
   \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_{i+1/2,j} \cdot \hat{\mathbf{n}}_{ii,kw} \: A_{i+1/2,j} & = &
   [\rho D_\alpha]_{i+1/2,j} \frac{Y_{\alpha i+1,j}-Y_{ \alpha i,j}}{\Delta x_{f(i)}} \Delta y_{c(j)} \label{eq:diffjj}
\end{eqnarray}
%
of opposite sign and same magnitude, as required by conservation.


\subsubsection*{Constant density diffusive flux for implicit time integration}

Now consider $ \mathbf{J_{d \alpha}'} = - D_\alpha \boldsymbol{\nabla} \left(\rho Y_\alpha \right)$, used in implicit integration of species transport. For cell $ii$, the first right hand side term of equation~\eqref{eq:intconvdiff} is
%
\begin{equation}
  \int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot \left(  \mathbf{J_{d \alpha}'} \right)  } d \Omega =
  \int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot \left( -  D_\alpha \boldsymbol{\nabla} \left( \rho Y_\alpha \right) \right)  } d \Omega =
  \int_{\partial \Omega_{ii}} {\left( -  D_\alpha \boldsymbol{\nabla} \left( \rho Y_\alpha \right) \right) \cdot \hat{\mathbf{n}}_{ii} } \: d \partial \Omega
\end{equation}
%
and in discrete form, we have
%
\begin{equation}
    \int_{\partial \Omega_{ii}} {\left( -  D_\alpha \boldsymbol{\nabla} \left( \rho Y_\alpha \right) \right) \cdot \hat{\mathbf{n}}_{ii} } \: d \partial \Omega = \sum^{nf_c}_{k=1}
    \left( -  D_\alpha \boldsymbol{\nabla} \left( \rho Y_\alpha \right) \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k \label{eq:discfvdiff2}
\end{equation}
%
where the last term in parenthesis is the average value of $ \mathbf{J_{d \alpha}'} $ over face $k$. Following similar derivation as in the previous section, we have for the east face ($ii,k_e \rightarrow i+1/2,j$) of a generic 2D cell:
%
\begin{eqnarray}
   \left( - D_\alpha \boldsymbol{\nabla} \left( \rho Y_\alpha \right) \right)_{i+1/2,j} \cdot \hat{\mathbf{n}}_{ii,ke} \: A_{i+1/2,j} & = &
   -[D_\alpha]_{i+1/2,j} \frac{ \left( \rho Y_\alpha \right)_{i+1,j} - \left( \rho Y_\alpha \right)_{i,j} } {\Delta x_{f(i)}} \Delta y_{c(j)}  \nonumber \\
   & = &  K_x \:\left( \rho Y_\alpha \right)_{i,j}-K_x \: \left( \rho Y_\alpha \right)_{i+1,j} \label{eq:diffii2}
\end{eqnarray}
where $K_x=[\D_\alpha]_{i+1/2,j} \frac{\Delta y_{c(j)}}{\Delta x_{f(i)}}$. By looking at equation~\eqref{eq:discfvdiff2}, we note that this term will add to the row $ii$ of a global diffusion matrix $\left[ \mathbf{A}_{diff \alpha} \right]$ corresponding to cell $(i,j)$ with single index $ii$. Now, for the adjacent cell indexed by $(i+1,j)$ with unique index $jj$, we have a term
%
\begin{equation}
   - \left[ -[D_\alpha]_{i+1/2,j} \frac{\left( \rho Y_\alpha \right)_{i+1,j} - \left( \rho Y_\alpha \right)_{i,j}}{\Delta x_{f(i)}} \Delta y_{c(j)} \right] = -K_x \: \left( \rho Y_\alpha \right)_{i,j}+K_x \: \left( \rho Y_\alpha \right)_{i+1,j} \label{eq:diffjj2}
\end{equation}
%
which will be added to the $jj$ equation of matrix $\left[ \mathbf{A}_{diff \alpha} \right]$, corresponding to the diffusive flux computed on face $i+1/2,j$. Therefore, the contribution of $x$-faces to $\left[ \mathbf{A}_{diff \alpha} \right]$ can be computed by assembling into the global matrix a set of face $2 \times 2$ matrices
%
\begin{equation}
   K^x_{i+1/2,j} = \left[ \begin{array}{cc} K_x & -K_x \\ -K_x & K_x \end{array} \right]
\end{equation}
whose coefficients map from local $1,2$ indexes to global $ii,jj$ indexes.
For the north ($i,j+1/2$) $y$-face the corresponding local diffusion matrix is
%
\begin{equation}
   K^y_{i,j+1/2} = \left[ \begin{array}{cc} K_y & -K_y \\ -K_y & K_y \end{array} \right]
\end{equation}
%
where $K_y=[D_\alpha]_{i,j+1/2} \frac{\Delta x_{c(i)}}{\Delta y_{f(j)}}$. Then, the global $\left[ \mathbf{A}_{diff \alpha} \right]$ is assembled looping through faces in the $x$ and $y$ directions. The same process would give a local matrix in the $z$ direction for 3D problems. The special case of domain boundary faces is discussed at the end of this section. The treatment of cell faces of irregular cut-cells that belong to the immersed geometry boundaries will be discussed in Section~\ref{sec:cc}.

\subsubsection*{FV spatial discretization of Advective Term}

The advective term of equation~\eqref{eq:intconvdiff} for cell $ii$ is discretized as
%
\begin{equation}
  \int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot  \left(  \rho Y_\alpha \mathbf{u} \right) } d \Omega =
  \int_{\partial \Omega_{ii}} { \left( \rho Y_\alpha \mathbf{u} \right) \cdot \hat{\mathbf{n}}_{ii} } \: d \partial \Omega =
  \sum^{nf_c}_{k=1} \left( \rho Y_\alpha \mathbf{u} \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k \label{eq:fvadv}
\end{equation}
%
where the velocity vector $\mathbf{u}$ can either refer to the advection velocity $\mathbf{v}$ or $\mathbf{v}'$. For the east ($ii,k_e \rightarrow i+1/2,j$)  face of a generic 2D cell:
%
\begin{eqnarray}
  \left( \rho Y_\alpha \mathbf{u} \right)_{i+1/2,j} \cdot \hat{\mathbf{n}}_{ii,ke} \: A_{i+1/2,j} & = &
  \overline{\left[ \rho Y_\alpha \right]}_{i+1/2,j} u_{i,j} \: \Delta y_{c(j)} \nonumber \\
  & = & \left(a^-_{i+1/2,j} \rho_{i,j} Y_{\alpha i,j} + a^+_{i+1/2,j} \rho_{i+1,j} Y_{\alpha i+1,j} \right) u_{i,j} \: \Delta y_{c(j)} \nonumber \\
  & = & C^-_x  \left( \rho Y_{\alpha} \right)_{i,j}  +  C^+_x  \left( \rho Y_{\alpha} \right)_{i+1,j} \label{eq:advfvcx}
\end{eqnarray}
%
where
%
\begin{eqnarray}
   C^-_x & = &  a^-_{i+1/2,j} \: u_{i,j}  \: \Delta y_{c(j)} \\
   C^+_x & = &  a^+_{i+1/2,j} \: u_{i,j} \:  \Delta y_{c(j)}
\end{eqnarray}
%
The coefficients $a^-_{i+1/2,j}$ and $a^+_{i+1/2,j}$ depend on the flux limiter chosen in the following form:
%
\begin{eqnarray}
   a^-_{i+1/2,j}  & = & \frac{1}{2} \left[ 1 + sgn(u_{i,j}) \left(1-B(r) \right) \right] \\
   a^+_{i+1/2,j} & = & \frac{1}{2} \left[ 1 - sgn(u_{i,j}) \left(1-B(r) \right) \right]
\end{eqnarray}
%
where $sgn(u_{i,j})$ is the sign of the face velocity and $B(r)=0$ implies \texttt{GODUNOV\_LIMITER} flux limiter, and $B(r)=1$ implies linear interpolation in cut-cell region faces (i.e. \texttt{CENTRAL\_LIMITER}).

Expression~\eqref{eq:advfvcx} defines the coefficients that will be added to row $ii$ of the global matrix $\left[ \mathbf{A}_{adv \alpha} \right]$. For cell $i+1,j$ indexed by single index $jj$ the corresponding term is
%
\begin{eqnarray}
  \left( \rho Y_\alpha \mathbf{u} \right)_{i+1/2,j} \cdot \hat{\mathbf{n}}_{jj,ke} \: A_{i+1/2,j} & = & - \overline{\left[ \rho Y_\alpha \right]}_{i+1/2,j} u_{i,j} \: \Delta y_{c(j)} \nonumber \\
  & = & -C^-_x \left( \rho Y_{\alpha} \right)_{i,j} - C^+_x \left( \rho Y_{\alpha} \right)_{i+1,j} \label{eq:advfvcxjj}
\end{eqnarray}
%
Leading to the face ${i+1/2,j}$ contribution to the advective matrix $\left[ \mathbf{A}_{adv \alpha} \right]$
%
\begin{equation}
   C^x_{i+1/2,j} = \left[ \begin{array}{cc} C^-_x & C^+_x \\ -C^-_x & -C^+_x \end{array} \right]
\end{equation}
%
In similar manner for the north ($i,j+1/2$) $y$-face, the corresponding local advection matrix is
%
\begin{equation}
   C^y_{i,j+1/2} = \left[ \begin{array}{cc} C^-_y & C^+_y \\ -C^-_y & -C^+_y \end{array} \right]
\end{equation}
%
where $C^-_y = a^-_{i,j+1/2} \: v_{i,j}  \: \Delta x_{c(i)}$, $C^+_y = a^+_{i,j+1/2} \: v_{i,j}  \: \Delta x_{c(i)}$.


\subsubsection*{FV treatment of boundary conditions} \label{sec:FVBCsCart}

%We look how boundary conditions affect each term in the discrete equation for boundary cells. We will maintain the guard cell approach (centered interpolation and difference on the boundary face), although using the prescribed face value or face derivative (first order one sided difference or face value) leads to the same result.

We are interested in \texttt{SOLID\_BOUNDARY} conditions. This boundary condition specifies zero mass transfer at the solid boundary. Kinematically, is specifies zero normal velocity at the solid boundary. Therefore the total mass flux integral for species $\alpha$ on a face $w$ belonging to a solid boundary adjacent to cell $g$ is (see figure~\ref{Fig:FVdisc}b):
%
\begin{equation}
    \dot{m}''_{\alpha,w} \; A_w = \left( \rho Y_\alpha \mathbf{v} - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right) \cdot \hat{\mathbf{n}}_{g,w} \: A_w = 0 \label{eq:fvsolidbc}
\end{equation}
%
where
\begin{equation}
    \dot{m}''_{\alpha,w} = \left( \rho Y_\alpha v_n - \rho D_\alpha \frac{\partial Y_\alpha}{\partial x_n} \right) = 0 \label{eq:fvsolidbc2}
\end{equation}
and from kinematics the normal advection velocity $\mathbf{v} \cdot \hat{\mathbf{n}}_{g,w}=v_n=0$, and from zero diffusive mass flux $\frac{\partial Y_\alpha}{\partial x_n}=0$. We note this last homogeneous Neumann condition on $Y_\alpha$ within our discretization as
%
\begin{equation}
   Y_{\alpha,w} = Y_{\alpha,g}
\end{equation}
%
The no mass flux condition at solid boundaries essentially specifies a zero component added to $\left\{ \mathbf{F_{BC \alpha}} \right\}$ on our previous equations, and no discretization matrix entries due to boundary conditions in case of implicit integration.

When using immersed boundary forcing to specify the no slip boundary conditions, a small value of normal velocity $v_{n,IB}$ (positive outside of cell $g$) is expected at face $w$. This is a consequence of the method used to reconstruct velocities at the boundary and the projection method used for time integration. Assuming a small known $v_{n,IB}$ and zero diffusive flux at face $w$, a transpiration mass flux:
%
\begin{equation}
    \rho_w Y_{\alpha,w} v_{n,IB}  \: A_w
\end{equation}
%
is added to $\left\{ \mathbf{F_{BC \alpha}} \right\}$ with corresponding sign.
%Note that this essentially explicit approximation of transpiration fluxes is also taken in implicit integration.


% Here, once implemented, also add SPECIFIED_MASS_FLUX BC.



\subsubsection{Cut-Cell Approaches} \label{sec:cc}

The treatment of irregular gas phase Eulerian grid faces and cells that remain after the inclusion of immersed bodies is described in this section.
We assume (see figure~\ref{Fig:FVdiscCC}a) the cut-cell definition algorithm has been successful in defining the following:
%
\begin{itemize}
   \item All \texttt{GASPHASE} faces in the $x$ and $y$ directions, which are tagged as regular or cut-faces.
           For these, we also know vertex points, areas (length in 2D), and face centroid location. We know for faces that are \texttt{GASPHASE}
           cut-faces the $i,j$ coordinates of the face they belong to.
  \item All \texttt{INBOUNDARY} faces, which arise from the intersection of the immersed bodies surface elements and the Eulerian
           grid cells. For these we know vertex points, areas (length in 2D), and face centroid location.
   \item All \texttt{GASPHASE} cut cells and regular cells. We know which faces define their boundary. For cut-cells, a list of cut-faces
           (type \texttt{GASPHASE} and \texttt{INBOUNDARY}) is provided. Also we know their volume and face centroid location.
  %\item By boolean difference of the Eulerian grid set and \texttt{GASPHASE} cut/regular cells and faces, we know the ~\texttt{SOLID} cells and faces on the grid (\textit{we might be interested on the \texttt{SOLID} cells and faces}).
\end{itemize}
%
%
\begin{figure}[h]
      %\centering
      \includegraphics[trim = 65mm 50mm 70mm 40mm, clip,
       width=0.49\linewidth]{../../../fig/fds/Geometry_Figures/CutCellsSketch.png}
      \includegraphics[trim = 65mm 40mm 50mm 30mm, clip,
       width=0.49\linewidth]{../../../fig/fds/Geometry_Figures/CutCellsBCSketch.png}
      \put(-350,-10){(a)}
      \put(-120,-10){(b)}
      \caption{Cut-cell FV discretization in 2D: (a) \texttt{GASPHASE} Cut-cell  $ii$ is bounded by: $\mathbf{1,2}$ \texttt{INBOUNDARY} faces, $\mathbf{3,6}$ \texttt{GASPHASE} cut-faces, and $\mathbf{4,5}$  \texttt{GASPHASE} regular faces.  (b) Geometrical elements used in the discretization of the diffusive term for \texttt{GASPHASE} cut-face $\mathbf{3}$, and \texttt{INBOUNDARY} face $\mathbf{1}$.}
	\label{Fig:FVdiscCC}
\end{figure}
%
\subsubsection{CC discretization of Diffusive term}  \label{sec:CCdiff}

%Here we discuss:
%\begin{itemize}
%  \item How are \texttt{GASPHASE} cut faces treated.
%  \item How are \texttt{INBOUNDARY} faces treated. Changes in diffusion matrix and $\mathbf{F_{bc}}$ for non homogeneous Neuman and Dirichlet.
%\end{itemize}

Consider the FV discretization of the diffusive term (equation~\eqref{eq:discfvdiff}) on cut-cell $ii$ in figure~\ref{Fig:FVdiscCC}b. We have
%
\begin{equation}
    \int_{\partial \Omega_{ii}} { \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right) \cdot \hat{\mathbf{n}}_{ii} } \: d \partial \Omega = \sum^{nf_c=6}_{k=1}
    \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k \label{eq:discfvdiffcc}
\end{equation}
%

In case of the diffusive like flux definition used in implicit time integration the expression is:
%
\begin{equation}
    \int_{\partial \Omega_{ii}} { \left( - D_\alpha \boldsymbol{\nabla} \left(  \rho Y_\alpha \right) \right) \cdot \hat{\mathbf{n}}_{ii} } \: d \partial \Omega = \sum^{nf_c=6}_{k=1}
    \left( - D_\alpha \boldsymbol{\nabla} \left(  \rho Y_\alpha \right) \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k \label{eq:discfvdiffcc2}
\end{equation}
%
The $nf_c=6$ faces that compose the boundary of the cut-cell can be divided in:

\subsubsection*{A. Faces $k=\mathbf{4},\mathbf{5}$ are \textit{regular} \texttt{GASPHASE} faces:}
The treatment of these is as described in the previous sections, with the caveat that the diffusive flux computation now involves the cut-cell centroid location (i.e. to compute spatial derivatives of $Y_\alpha$ and interpolation of $\rho D_\alpha$ in equation~\eqref{eq:discfvdiffcc}).

\subsubsection*{B. Faces $k=\mathbf{3},\mathbf{6}$ are \texttt{GASPHASE} \textit{cut-faces}:}

In order to compute the discrete term of equation~\eqref{eq:discfvdiffcc} on these faces, the factor $\rho D_\alpha$ needs to be interpolated from cell centroids to the face centroid and the spatial derivative of $Y_\alpha$ computed at the face centroid. The use of face centroids is required to maintain spatial accuracy when going from the left to right hand side of~\eqref{eq:discfvdiffcc} in difference FV methods. To maintain accuracy of the overall cut cell method, interpolation and differentiation to the cut-face centroid must also be sufficiently accurate.

The following discussion is fairly general, in the sense that it is agnostic to the interpolation methods employed (Lagrange polynomials, isoparametric, least squares, etc.).
Consider a stencil of points $e=1,...,ne$, where the scalar $Y_\alpha$ is assumed defined (i.e. cell centroids). Then for a location of interest $ck$ (i.e. the centroid of \texttt{GASPHASE} cut-face $k=3,6$), interpolation and derivatives of $Y_\alpha$ at $ck$ can be obtained as
%
\begin{eqnarray}
   Y_\alpha(\mathbf{x}_{ck}) & = & \sum^{ne}_{e=1} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) Y_{\alpha e} \label{eq:interpqfluid} \\
   \frac{\partial Y_\alpha(\mathbf{x}_{ck})}{\partial x_i} & = & \sum^{ne}_{e=1} \frac{\partial \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e)}{\partial x_i} Y_{\alpha e} \label{eq:interpdqfluid}
\end{eqnarray}
%
where $\mathbf{x}_e$ is the location of stencil point $e$, $Y_{\alpha e} =Y_\alpha(\mathbf{x}_e)$  and $\phi_e(\mathbf{x}-\mathbf{x}_e)$, $e=1,...,ne$ is a suitable set of interpolation functions.
For this discussion we assume that no boundary values of $Y_\alpha$ are involved on the interpolation. For these faces the diffusive flux is:
%
\begin{eqnarray}
  \left( - \rho D_q \boldsymbol{\nabla} Y_\alpha \right)_k & = & - \sum^{ne}_{m=1} \phi_m(\mathbf{x}_{ck}-\mathbf{x}_m) [\rho D_\alpha]_m \times
      \sum^{ne}_{e=1} \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \; Y_{\alpha e} \nonumber \\
      &=& - \sum^{ne}_{e=1} [\rho D_\alpha]_k \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \; Y_{\alpha e}
\end{eqnarray}
%
where
%
\begin{equation}
 [\rho D_\alpha]_k = \left[ \sum^{ne}_{m=1} \phi_m(\mathbf{x}_{ck}-\mathbf{x}_m) [\rho D_\alpha]_m \right]
\end{equation}
%
The contribution to the equation of the cut-cell $jj$ on the \textit{low} side of the  \texttt{GASPHASE} cut-cell is, following what was seen for regular cells:
%
\begin{equation}
   \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{jj,k} \: A_k = - \sum^{ne}_{e=1} K_{e,k} \; Y_{\alpha e} \label{eq:kekdiff}
\end{equation}
%
where $K_{ek} = [\rho D_\alpha]_k \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \cdot \hat{\mathbf{n}}_{jj,k} \: A_k$.

Similarly for constant density diffusive fluxes (equation~\eqref{eq:discfvdiffcc2}) on a \texttt{GASPHASE} cut-face $k$:
%
\begin{equation}
   \left( - D_\alpha \boldsymbol{\nabla} \left( \rho Y_\alpha \right) \right)_k \cdot \hat{\mathbf{n}}_{jj,k} \: A_k = - \sum^{ne}_{e=1} K_{e,k}' \; \left( \rho Y_{\alpha} \right)_e \label{eq:kekdiff2}
\end{equation}
%
where $K_{ek}' = [D_\alpha]_k \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \cdot \hat{\mathbf{n}}_{jj,k} \: A_k$.

We see by this last expression that unknown connectivity on the row $jj$ of $\left[ \mathbf{A}_{diff \alpha} \right]$ has been increased, that is, now it involves $ne$ unknowns instead of 2. For cut-cell $ii$ on the high side of cut-face $k$ the contribution is $\sum^{ne}_{e=1} K_{e,k} \; Y_{\alpha e}$ or  $\sum^{ne}_{e=1} K_{e,k}' \; \left( \rho Y_{\alpha} \right)_e$ (different sign).


%%
%%
\subsubsection*{1D Linear Isoparametric Polynomial Interpolation:}
%%%%%

Consider a stencil with cells $jj$ and $ii$ as in figure~\ref{Fig:FVdiscCC}b. Their centroids are located at $\mathbf{x}_{jj}$ and $\mathbf{x}_{ii}$ respectively. Parameterizing point locations along the segment that unites $\mathbf{x}_{jj}$ to $\mathbf{x}_{ii}$ by $-1 \le \xi \le 1$, we have
%
\begin{eqnarray}
   \mathbf{x}(\xi)  & = &  \phi_1(\xi) \mathbf{x}_{jj} + \phi_2(\xi) \mathbf{x}_{ii} \label{eq:xxi} \\
   Y_\alpha (\xi)  & = & \phi_1(\xi) Y_{\alpha jj} + \phi_2(\xi) Y_{\alpha ii}
\end{eqnarray}
%
\begin{equation}
[ \rho D_\alpha ] (\xi)  =  \phi_1(\xi) [ \rho D_\alpha ]_{jj} + \phi_2(\xi) [ \rho D_\alpha ]_{ii}
\end{equation}
where
%
\begin{equation}
   \phi_1(\xi) = \frac{1}{2} \left(1 - \xi \right) \; ; \; \phi_2(\xi) = \frac{1}{2} \left(1 + \xi \right)
\end{equation}
%
We assume the value of the parameter where the cut-face lies known $\xi_p$, corresponding to the face and centroid segment intersection $\mathbf{x}_p$. From equation~\eqref{eq:xxi} with $\mathbf{x}=x\hat{\mathbf{i}}+y\hat{\mathbf{j}}$, we have
%
\begin{eqnarray}
   \frac{\partial \xi}{\partial x} = \frac{2}{(x_{ii}-x_{jj})} \; ; \;  \frac{\partial \xi}{\partial y} = \frac{2}{(y_{ii}-y_{jj})}
\end{eqnarray}
%
and the interpolation function gradient components are:
%
\begin{eqnarray}
  \frac{\partial\phi_1}{\partial x} & = & \frac{\partial\phi_1}{\partial \xi} \frac{\partial \xi}{\partial x}=-\frac{1}{(x_{ii}-x_{jj})} \\
  \frac{\partial\phi_1}{\partial y} & = & \frac{\partial\phi_1}{\partial \xi} \frac{\partial \xi}{\partial y}=-\frac{1}{(y_{ii}-y_{jj})} \\
  \frac{\partial\phi_2}{\partial x} & = & \frac{\partial\phi_2}{\partial \xi} \frac{\partial \xi}{\partial x}= \frac{1}{(x_{ii}-x_{jj})} \\
  \frac{\partial\phi_2}{\partial y} & = & \frac{\partial\phi_2}{\partial \xi} \frac{\partial \xi}{\partial y}= \frac{1}{(y_{ii}-y_{jj})}
\end{eqnarray}
%
which are constant. Finally, being $k=3$, for cut-cell $jj$ the normal is $\hat{\mathbf{n}}_{jj,k}=\hat{\mathbf{j}}$, and expression~\eqref{eq:kekdiff} reduces to
%
\begin{equation}
  \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{jj,k} \: A_k =
  -\left[  -K_{jj,k} \: Y_{\alpha jj} +  K_{ii,k} \: Y_{\alpha ii} \right]
\end{equation}
%
where $K_{jj,k}=K_{ii,k}=[ \rho D_\alpha ]_p \frac{A_k}{(y_{ii}-y_{jj})}$. Note that, an assumption of approximation of the cut-face centroid location by point $p$ has been made. For cut-cell $ii$ on the high side of face $k=3$ the expression is
%
\begin{equation}
  \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k =
  \left[  -K_{jj,k} \: Y_{\alpha jj} +  K_{ii,k} \: Y_{\alpha ii} \right]
\end{equation}
%

Similar expressions would be found for a \texttt{GASPHASE} cut-face normal to the $x$ direction. The last two expressions are the cut-cell versions of what was seen on Section~\ref{Sec:FVdiff}. Similar expressions are found for constant density diffusive fluxes, where $K_{jj,k}'=K_{ii,k}'=[D_\alpha ]_p \frac{A_k}{(y_{ii}-y_{jj})}$ coefficients are associated with variables $\left( \rho Y_{\alpha} \right)_{ii}$, and $\left( \rho Y_{\alpha} \right)_{jj}$.


\subsubsection*{C. Faces $\mathbf{1},\mathbf{2}$ are \texttt{INBOUNDARY} \textit{cut-faces}:}

Consider \texttt{INBOUNDARY} cut-face $k=\mathbf{1}$ in figure~\ref{Fig:FVdiscCC}b. If a Neumann immersed boundary condition is specified in such a face, the term
%
\begin{equation}
  \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k = [\rho D_\alpha]_k \partial \gamma_{ii,k} \: A_k
\end{equation}
%
where as before, the prescribed normal derivative of $Y_\alpha$ on this face $\partial Y_\alpha / \partial n |_k=\partial  \gamma_{ii,k}$ is defined in the direction of $\hat{\mathbf{n}}_1$.
This term adds to $\mathbf{F_{BC \alpha}}$ in the location of cell $ii$. Similarly to what was seen on section~\ref{sec:FVBCsCart} for \texttt{SOLID\_BOUNDARIES}, zero diffusive mass flux implies a homogeneous $\partial Y_\alpha / \partial n |_k=0$ for the mass fraction field of species $\alpha$.



% We note that, as seen for domain boundary faces, there is no modification of the diffusion matrix by \texttt{INBOUNDARY} cut-faces.

%When using \textit{Dirichlet} boundary conditions for the scalar $Y_\alpha$ on immersed geometries, an estimation of $\partial Y_\alpha / \partial n |_k$ is required. The data is $Y_{\alpha k}=\gamma_{ii,k}$, and we assume we know its values on the face vertices (i.e. $Y_{\alpha lk}$ in $\mathbf{x}_{lk}$, and $Y_{\alpha hk}$ in $\mathbf{x}_{hk}$ where here $k=1$). Therefore, the $\partial Y_\alpha / \partial n |_k$ is generally approximated using a stencil composed of: (a) a set of $ne$ points corresponding to \texttt{GASPHASE} cut-cell centroids, and (b) a set of $nik$ points located on the \texttt{INBOUNDARY} cut-face $k$. The corresponding expressions for interpolation and derivative approximations are
%%
%\begin{eqnarray}
%   Y_\alpha(\mathbf{x}) & = &  \sum^{ne}_{e=1} \phi_e(\mathbf{x}-\mathbf{x}_e) Y_{\alpha e}
%                              +   \sum^{ne+nik}_{e=ne+1} \phi_e(\mathbf{x}-\mathbf{x}_e) \gamma_{ek}  \label{eq:interpqfluidbdy} \\
%   \frac{\partial Y_\alpha(\mathbf{x})}{\partial x_i} & = &  \sum^{ne}_{e=1} \frac{\partial \phi_e(\mathbf{x}-\mathbf{x}_e)}{\partial x_i} Y_{\alpha e}
%                                                                        + \sum^{ne+nik}_{e=ne+1} \frac{\partial \phi_e(\mathbf{x}-\mathbf{x}_e)}{\partial x_i} \gamma_{ek}   \label{eq:interpdqfluidbdy}
%\end{eqnarray}
%%
%where $\gamma_{ek}$, $e=1,...,nik$, are the Dirichlet data on the \texttt{INBOUNDARY} cut-face. Assuming derivatives approximation to the cut-cell centroid $ck$, the diffusive flux on the face is
%%
%\begin{eqnarray}
%  \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k =  - [\rho D_\alpha]_k \sum^{ne}_{e=1} \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) Y_{\alpha e}
%                                                                                - [\rho D_\alpha]_k \sum^{ne+nik}_{e=ne+1} \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e)\gamma_{ek}
%\end{eqnarray}
%%
%and
%\begin{eqnarray}
%  \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k =
%  &-&  \sum^{ne}_{e=1} [\rho D_\alpha]_k \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \cdot \hat{\mathbf{n}}_{ii,k} \:  A_k  \: Y_{\alpha e}  \nonumber \\
%  &-&[\rho D_\alpha]_k \left[ \sum^{ne+nik}_{e=ne+1} \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \gamma_{ek} \right] \cdot \hat{\mathbf{n}}_{ii,k} \:  A_k
%\end{eqnarray}
%%
%But as $\hat{\mathbf{n}}_{ii,k}=-\hat{\mathbf{n}}_{k}$,
%%
%\begin{equation}
%  \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k = \sum^{ne}_{e=1} K_{e,k} \: Y_{\alpha e} + F_{k}
%\end{equation}
%%
%with $F_{k}=[\rho D_\alpha]_k \left[ \sum^{ne+nik}_{e=ne+1} \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \gamma_{ek} \right] \cdot \hat{\mathbf{n}}_{k} \:  A_k$ being added to $\mathbf{F_{bc}}$, and $K_{e,k}=[\rho D_\alpha]_k \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \cdot \hat{\mathbf{n}}_{k} \:  A_k$ are coefficients to be assembled into row $ii$ of $\left[ \mathbf{A}_{diff \alpha} \right]$. As found for domain boundary faces, Dirichlet conditions on \texttt{INBOUNDARY} cut-faces modify both diffusion matrix and boundary condition vector.
%
%
%
%
%
%%%
%%%
%\subsubsection*{2D Linear Isoparametric Polynomial Interpolation:}
%%%%%%
%In figure~\ref{Fig:FVdiscCC}b, we define the triangle $T_{ii,k}$, connecting vertices $\mathbf{x}_{ii}$, $\mathbf{x}_{hk}$ and $\mathbf{x}_{lk}$. These vertices take local triangle indexes $1$, $2$, and $3$ respectively. Data values of $Y_\alpha$ on face $k=1$ are given ($Y_{\alpha hk}=\gamma_{hk}$ and $Y_{\alpha lk}=\gamma_{lk}$). The isoparamteric interpolation formulas on a unitary master triangle are
%%
%\begin{eqnarray}
%  \mathbf{x}(\xi,\eta) & = & \phi_1(\xi)  \mathbf{x}_{ii} + \phi_2(\eta) \mathbf{x}_{hk} + \phi_3(\xi,\eta) \mathbf{x}_{lk} \\
%  Y_\alpha(\xi,\eta) & = & \phi_1(\xi)  Y_{\alpha ii} + \phi_2(\eta) Y_{\alpha hk} + \phi_3(\xi,\eta)Y_{\alpha lk}
%\end{eqnarray}
%%
%where $\phi_1(\xi)=\xi$, $\phi_2(\eta)=\eta$, and $\phi_3(\xi,\eta)=1-\xi-\eta$. The gradient of a function respect to the spatial coordinates is related  to the derivatives respect to natural coordinates $\xi,\eta$ as $\boldsymbol{\nabla} () = [\mathbf{J}]^{-1} [\partial/\partial \xi \; \partial/\partial \eta]^T$. The Jacobian matrix and its inverse are
%%
%\begin{equation}
%  [\mathbf{J}] = \left[ \begin{array}{cc} x_1-x_3 & y_1-y_3 \\ x_2-x_3 & y_2-y_3 \end{array} \right] \; ; \;
%  [\mathbf{J}]^{-1} = \frac{1}{det([\mathbf{J}])} \left[ \begin{array}{cc} y_2-y_3 & y_3-y_1 \\ x_3-x_2 & x_1-x_3 \end{array} \right]
%\end{equation}
%%
%where $det([\mathbf{J}])$ is the determinant of $[\mathbf{J}]$. As $\phi_1$ is the interpolation function related to centroid $ii$ we have
%$\boldsymbol{\nabla} \phi_1 = det([\mathbf{J}])^{-1} [(y_2-y_3) \; (x_3-x_2)]^T$ (constant), and the unique factor that will add to the diagonal of $\left[ \mathbf{A}_{diff \alpha} \right]$ in row $ii$ is
%%
%\begin{equation}
%  K_{e,k}=[\rho D_\alpha]_k det([\mathbf{J}])^{-1} \left[ (y_2-y_3)n^x_k + (x_3-x_2) n^y_k \right] A_k
%\end{equation}
%%
%while the coefficient adding to $\mathbf{F^{diff}_{bc}}$ is
%%
%\begin{equation}
%  F_{k}=[\rho D_\alpha]_k \left[ \boldsymbol{\nabla} \phi_2 \gamma_{hk} + \boldsymbol{\nabla} \phi_3 \gamma_{lk} \right] \cdot \hat{\mathbf{n}}_{k} \:  A_k
%\end{equation}
%%
%where $\hat{\mathbf{n}}_{k}=n^x_k \hat{\mathbf{i}}+n^y_k\hat{\mathbf{j}}$, and
%%
%\begin{eqnarray}
%   \boldsymbol{\nabla} \phi_2 & = & det([\mathbf{J}])^{-1} [(y_3-y_1) \; (x_1-x_3)]^T \\
%   \boldsymbol{\nabla} \phi_3 & = & -(\boldsymbol{\nabla} \phi_1+\boldsymbol{\nabla} \phi_2)
%\end{eqnarray}

\subsubsection{CC discretization of Advective term}

We study the discretization of the advective term~\eqref{eq:fvadv} for the six faces of cell $ii$ in figure~\ref{Fig:FVdiscCC}b, making use of the interpolation and derivative approximation expressions~\eqref{eq:interpqfluid}-\eqref{eq:interpdqfluid}, seen on the previous section.

\subsubsection*{A. Faces $k=\mathbf{4},\mathbf{5}$ are \textit{regular} \texttt{GASPHASE} faces:}
The treatment of these is as described in the previous section, with the caveat that the advective flux computation for $\rho Y_\alpha$ may now involve the cut-cell centroid location in the case of spatial interpolation to the face (i.e. \texttt{CENTRAL\_LIMITER} flux limiter).

\subsubsection*{B. Faces $k=\mathbf{3},\mathbf{6}$ are \texttt{GASPHASE} \textit{cut-faces}:}

For \texttt{GASPHASE} cut-face $k=3$ in figure~\ref{Fig:FVdiscCC}b, the advective term corresponding to cell $jj$ is
%
\begin{equation}
  \left( \rho Y_\alpha \mathbf{u} \right)_k \cdot \hat{\mathbf{n}}_{jj,k} \: A_k = \overline{[\rho Y_\alpha]}_k \left( \mathbf{u}_k \cdot \hat{\mathbf{n}}_{jj,k} \right) \: A_k \label{eq:convgcutface}
\end{equation}
%
where the average face normal velocity $\mathbf{u}_k \cdot \hat{\mathbf{n}}_{jj,k} $ (i.e. approximated by the face centroid value) is assumed known, and $\overline{[\rho Y_\alpha]}_k$ is obtained by flux limited interpolation to the face.

Consider a flux limited spatial interpolation of the form~\eqref{eq:interpqfluid} for $[\rho Y_\alpha]$ to face centroid $ck$
%
\begin{equation}
  \overline{[\rho Y_\alpha]}_k = \sum^{ne}_{e=1} \overline{\phi}_e(\mathbf{x}_{ck}-\mathbf{x}_e) \left( \rho Y_{\alpha} \right)_e
\end{equation}
%
where $e=1,...,ne$ refers to a suitable flux limited stencil of \texttt{GASPHASE} cell centroids. The $\overline{\phi}_e$ are the flux limited interpolation functions  related to said centroids. Then inserting this approximation in equation~\eqref{eq:convgcutface}
%
\begin{equation}
  \left( \rho Y_\alpha \mathbf{u} \right)_k \cdot \hat{\mathbf{n}}_{jj,k} \: A_k = \sum^{ne}_{e=1}  C_{e,k} \left( \rho Y_{\alpha} \right)_e
\end{equation}
being $C_{e,k}=\overline{\phi}_e(\mathbf{x}_{ck}-\mathbf{x}_e) \left( \mathbf{u}_k \cdot \hat{\mathbf{n}}_{jj,k} \right) \: A_k$.


%%
%%
\subsubsection*{1D Linear Isoparametric Polynomial Interpolation:}
%%%%%

For face $k=3$ consider the $y$ coordinate interpolation to point $p$ using cell centroids $\mathbf{x}_{jj}$ and $\mathbf{x}_{ii}$ (see  figure~\ref{Fig:FVdiscCC}b)
%
\begin{equation}
  y_p =  \phi_1(\xi_p) y_{jj} + \phi_2(\xi_p) y_{ii}
\end{equation}
%
As we know $y_p$ (face coordinate), the previous expression is used to obtain the natural coordinate directly
%
\begin{equation}
  \xi_p = \frac{y_p - \frac{1}{2} \left( y_{jj} + y_{ii} \right)}{\frac{1}{2} \left( y_{ii} - y_{jj} \right)}
\end{equation}
%
we see that if $y_p$ is located in the midpoint between $y_{jj}$ and $y_{ii}$, then $\xi_p=0$, the interpolation functions $\phi_1(\xi_p)=\phi_2(\xi_p)=1/2$, and the \textit{central} interpolation is recovered. In similar manner as seen for diffusion, for cell $jj$
$\hat{\mathbf{n}}_{jj,k}=\hat{\mathbf{j}}$ and:
%
\begin{equation}
   \left( \rho Y_\alpha \mathbf{u} \right)_k \cdot \hat{\mathbf{n}}_{jj,k} \: A_k = C_{1,k} \left( \rho Y_\alpha \right)_{jj} + C_{2,k} \left( \rho Y_\alpha \right)_{ii}
\end{equation}
%
where $C_{1,k}=\phi_1(\xi_p) v_k \: A_k$, $C_{2,k}=\phi_2(\xi_p) v_k \: A_k$. For cell $ii$ the advective flux term is
($\hat{\mathbf{n}}_{jj,k}=-\hat{\mathbf{j}}$).
%
\begin{equation}
   \left( \rho Y_\alpha \mathbf{u} \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k = - C_{1,k}  \left( \rho Y_\alpha \right)_{jj} - C_{2,k}  \left( \rho Y_\alpha \right)_{ii}
\end{equation}
%
Similar expressions can be found for \texttt{GASPHASE} cut-faces aligned in other directions. For a \textit{Godunov} (1st order upwind) interpolation,
we replace $\phi_1(\xi_p)$ by $\overline{\phi_1}(v_k)=1/2*(1+sgn(v_k))$ and $\phi_2(\xi_p)$ by $\overline{\phi}_2(v_k)=1/2*(1-sgn(v_k))$ on the $C_{1,k}$ and $C_{2,k}$ definitions. Godunov flux limited interpolation is used in these cut-faces.


\subsubsection*{C. Faces $\mathbf{1},\mathbf{2}$ are \texttt{INBOUNDARY} \textit{cut-faces}:}

The treatment of advective terms in boundary cut-faces is the same as defined in section~\ref{sec:FVBCsCart}, in particular in the source arising from the immersed boundary treatment derived transpiration velocities.

\subsubsection{Unsteady evolution: Fully explicit time integration for scalars} \label{sec:exscl}

It is well known that cut-cell methods pose a significant time constraint when used with explicit time integration methods.
This is so because, inevitably for general problems there will arise \texttt{GASPHASE} cut-cells whose small-size will penalize
severely the time step. We recall that each cell on the \texttt{GASPHASE}, including cut-cells, needs to meet CFL and Von Neuman stability constraints. Several different ways have been proposed in the literature to deal with this problem, i.e. cell merging, mixing or linking methods. In general, these lead to ad hoc selection procedures for surrounding cells, having to deal with many special cases, and in some cases potential solution deterioration close to the boundary.

In our fully explicit time integration implementation we use a simple and robust procedure to address this problem. Within the scheme for numeration of scalar cell unknowns, a test is performed on cut-cells. If the cut-cell volume is less than the threshold volume $V_{thr}= c_{thr} V_{cart}$, where $V_{cart}$ is the local Cartesian cell volume and $c_{thr}<1$ is a threshold factor, the unknown number this cell takes is the one of an adjacent cell which has a volume larger than $V_{thr}$. This mathematically defines a single control volume of the two linked cells. Cell volumes are added in building the mass matrix for the FV discretization, and fluxes and matrix terms are added with their corresponding signs. Note that, flux quantities corresponding to the common face of these two cells effectively cancel on the single equation for the set.
Alternatively, if after a number of cell numbering iterations an unlinked small cell persists in the mesh, the method tries to link such cut-cell to the closest numbered regular cell. Finally, if this does not succeed the small cell is given the $V_{thr}$ volume. This last resort makes the scheme fault proof, albeit at an error penalty on these "trapped" small cells.



\subsubsection{Unsteady evolution: Explicit-implicit time integration for scalars} \label{sec:eximscl}


We define another way to handle the small cell-problem in a robust manner using implicit integration of scalars. Nevertheless, to reduce the cost of building discretization matrices and solving linear systems of the size of the vector $\{\mathbf{\rho Y_\alpha}\}$, that contains unknowns for the whole domain, we only solve implicitly our problem in a band of cells on the cut-cell region. See figure~\ref{Fig:DOMEXIM}. We gather the unknowns from  $\{\mathbf{\rho Y_\alpha}\}$ that belong to this implicit region in vector $\{\mathbf{\rho Y_\alpha}\}_{IM}$, and their complement on $\{\mathbf{\rho Y_\alpha}\}$ in the explicit region unknowns $\{\mathbf{\rho Y_\alpha}\}_{EX}$. In general, the explicit region is advanced first. For the integration to be conservative, the normal flux of $q$ used in the implicit-explicit boundary has to be the same for implicit and explicit region time advancement. This constraint provides the non-homogeneous flux boundary condition to be used in integrating the implicit region.

%
\begin{figure}[h]
      \centering
      \includegraphics[trim = 65mm 55mm 70mm 40mm, clip,
       width=0.49\linewidth]{../../../fig/fds/Geometry_Figures/DomainsEXIM.png}
      \caption{Explicit-implicit integration of \texttt{GASPHASE} region including cut-cells. The regular grid in $\Omega_{EX}$ is integrated explicitly, while a band region $\Omega_{IM}$ of three cells including cut-cells surrounding the solid is integrated implicitly. Fluxes computed explicitly in the boundary $\Gamma_{EXIM}$ are used as boundary condition for the implicit solve.}
	\label{Fig:DOMEXIM}
\end{figure}
%

The idea of using an explicit flux of $q$ definition to separate subdomains of integration is not new. In their landmark work applied to parabolic problems, Dawson and Dupont~\cite{Dawson:1994} used it to decompose the spatial domain in different subdomains that would be integrated implicitly. The computation speedup is obtained here when solving smaller matrix systems in parallel, but in general the stability of the method is reduced compared to solving the original implicit problem. This work has led to a large body of research and literature on what is called explicit-implicit domain decomposition (EIDD) or time-stepping. Recently, May and Berger~\cite{May:2014,May:2017} used similar arguments to couple a MUSCL solver on an explicit integration region with either implicit Euler or trapezoidal rule on the implicit region, with application to the linear advection equation. They prove strong stability preservation (SSP) for the MUSCL-implicit Euler combination, alleviating successfully the small cut-cell problem on their embedded boundary method. Although the SSP property is only required for convergence in nonlinear advection, in many linear advection cases is desirable to reduce numerical oscillation~\cite{Gottlieb:2001}.
In our case we are interested in explicit-implicit solve combinations for Runge-Kutta methods. A schematic of the method can be given considering the explicit-implicit Euler pair: given $\{\mathbf{\rho Y_\alpha}\}^n$, $\rho^n$, $\mathbf{v}^n$ and $D_\alpha^n$ known on $\Omega_{EX} \cup \Omega_{IM}$ at time instance $n$, and a suitable time step $\Delta t$
%
\begin{enumerate}
  \item Integrate the discrete version of equation~\eqref{eq:bal3} on $\Omega_{EX}$, $\{\mathbf{q}_{EX}\}$ using Explicit Euler:
  \begin{equation}
  \{\mathbf{\rho Y_\alpha}\}^{n+1}_{EX}=\{\mathbf{\rho Y_\alpha}\}^{n}_{EX} + \Delta t \mathbf{L}(\{\mathbf{\rho Y_\alpha}\}^n,t_n)
  \end{equation}
  Where $\mathbf{L}()$ corresponds to the discrete convection, diffusion and source terms.
  \item Using $\{\mathbf{\rho Y_\alpha}\}^n$, define fluxes $\left( [\rho Y_\alpha \mathbf{v} -\rho D_\alpha \boldsymbol{\nabla} Y_\alpha] \cdot \hat{\mathbf{n}} \right)^n$ at the $\Gamma_{EXIM}$ boundary.
  \item Using the previously defined scalar fluxes as boundary condition, perform integration of equation~\eqref{eq:bal4} on $\Omega_{IM}$ using implicit Euler:
  \begin{equation}
  \{\mathbf{\rho Y_\alpha}\}^{n+1}_{IM}=\{\mathbf{\rho Y_\alpha}\}^{n}_{IM}+\Delta t \mathbf{L}'(\{\mathbf{\rho Y_\alpha}\}^{n+1},t_{n+1})
  \end{equation}
\end{enumerate}
%
where $\mathbf{L}()'$ is the convection-diffusion operator plus source terms in equation~\eqref{eq:bal4}. We note that step 2, ensures the conservation of $\rho Y_\alpha$ on $\Omega_{EX} \cup \Omega_{IM}$.

We are interested in using as explicit scheme the strong stability preserving SSPRK2~\cite{Gottlieb:2001} method currently used in FDS.
Following~\cite{May:2014}, we use as implicit counterparts either implicit Euler or trapezoidal rule, the latter providing second order time accuracy, albeit at the loss of the SSP property for time step sizes corresponding to local CFL number greater than 2\textbf{?}~\cite{Ketcheson:2009,Gottlieb:2009}. As expected, this implicit scheme is still linearly stable for all time steps. Indeed, the time step restriction for the integration in both time stepping combinations is the restriction obtained for the explicit region.
We define the following schemes:

\subsubsection*{Explicit SSPRK2 + Implicit 1st order BE:}

%
\begin{enumerate}
   \item First Stage: Euler pair with time step $\Delta t$
   Forward Euler FE for $EX$ region and Backward Euler BE for $IM$ region. Boundary conditions for the implicit solve are provided by the explicit mass fluxes $\left( [\rho Y_\alpha \mathbf{v} -\rho D_\alpha \boldsymbol{\nabla} Y_\alpha] \cdot \hat{\mathbf{n}} \right)^n$ at the $\Gamma_{EXIM}$ boundary.
   \begin{eqnarray}
   \{\mathbf{\rho Y_\alpha}\}^*_{EX}&=&\{\mathbf{\rho Y_\alpha}\}^{n}_{EX} +
   \Delta t \mathbf{L}(\{\mathbf{\mathbf{\rho Y_\alpha}}\}^n,t_n) \label{eq:m1s1ex} \\
   \{\mathbf{\rho Y_\alpha}\}^*_{IM}&=&\{\mathbf{\rho Y_\alpha}\}^{n}_{IM} +
   \Delta t \mathbf{L}'(\{\mathbf{\mathbf{\rho Y_\alpha}}\}^{n+1},t_{n+1}) \label{eq:m1s1im}
   \end{eqnarray}
  \item Second Stage: SSPRK2 corrector for $EX$ and $IM$ regions. At the $IM$ region $L'()$ is evaluated implicitly. Boundary conditions for the implicit solve are given by $\left( [\rho Y_\alpha \mathbf{v} -\rho D_\alpha \boldsymbol{\nabla} Y_\alpha] \cdot \hat{\mathbf{n}} \right)^*$.
   \begin{eqnarray}
   \{\mathbf{\rho Y_\alpha}\}^{n+1}_{EX}&=& \frac{1}{2}\left( \{\mathbf{\rho Y_\alpha}\}^{n}_{EX} +
                                                                                              \{\mathbf{\rho Y_\alpha}\}^*_{EX} \right) +
                                                   \frac{\Delta t}{2} \mathbf{L}(\{\mathbf{\rho Y_\alpha}\}^*,t_{n+1})  \label{eq:m1s2ex} \\
   \{\mathbf{\rho Y_\alpha}\}^{n+1}_{IM}&=& \frac{1}{2}\left( \{\mathbf{\rho Y_\alpha}\}^{n}_{IM} +
                                                                                             \{\mathbf{\rho Y_\alpha}\}^{*}_{IM} \right) +
                                                   \frac{\Delta t}{2} \mathbf{L}'(\{\mathbf{\rho Y_\alpha}\}^{n+1},t_{n+1}) \label{eq:m1s2im}
   \end{eqnarray}
\end{enumerate}
%
In the hypothetical case where the $IM$ region covers the whole computational domain, it is easy to see that the implicit integrations on first and second stages correspond to the backward Euler method (replace~\eqref{eq:m1s1im} in~\eqref{eq:m1s2im}). This combination maintains the sequence of SSPRK2 used in FDS. It provides a first order accurate method which is also SSP for any time step, as both SSPRK2 and BE methods are, and the fluxes at the explicit-implicit boundary are evaluated as in SSPRK2.

\subsubsection*{Explicit SSPRK2 + Implicit 1st order BE and Trapezoidal rule:}
%
\begin{enumerate}
   \item First Stage: Same solution as in the previous method and also integrate explicitly the $IM$ region (in practice an explicit integration of the whole domain is done, followed by BE on the $IM$ region)
   \begin{equation}
     \{\mathbf{\rho Y_\alpha}\}^{ex}_{IM} = \{\mathbf{\rho Y_\alpha}\}^{n}_{IM}+\Delta t \mathbf{L}'(\{\mathbf{\rho Y_\alpha}\}^n,t_n) \label{eq:m2s1imex}
   \end{equation}

  \item Second Stage: SSPRK2 corrector for $EX$ and $IM$ regions. At the $IM$ region $\mathbf{L}'()$ is evaluated implicitly, boundary fluxes are defined as for the previous corrector, and $\{\mathbf{\rho Y_\alpha}\}^{*}_{IM}$ is replaced by $\{\mathbf{\rho Y_\alpha}\}^{ex}_{IM}$
   \begin{eqnarray}
   \{\mathbf{\rho Y_\alpha}\}^{n+1}_{EX}&=& \frac{1}{2}\left( \{\mathbf{\rho Y_\alpha}\}^{n}_{EX} +
                                                                                              \{\mathbf{\rho Y_\alpha}\}^*_{EX} \right) +
                                                   \frac{\Delta t}{2} \mathbf{L}(\{\mathbf{\rho Y_\alpha}\}^*,t_{n+1})  \label{eq:m2s2ex} \\
   \{\mathbf{\rho Y_\alpha}\}^{n+1}_{IM}&=& \frac{1}{2}\left( \{\mathbf{\rho Y_\alpha}\}^{n}_{IM} +
                                                                                             \{\mathbf{\rho Y_\alpha}\}^{ex}_{IM} \right) +
                                                   \frac{\Delta t}{2} \mathbf{L}'(\{\mathbf{\rho Y_\alpha}\}^{n+1},t_{n+1}) \label{eq:m2s2im}
   \end{eqnarray}
\end{enumerate}
%
Replacing equation~\eqref{eq:m2s1imex} in~\eqref{eq:m2s2im} we see that this last step correspond to the Trapezoidal rule applied to the $IM$ region.


\subsubsection{Implicit time discretization and simplifications}

%Consider the exact decomposition of the diffusive flux in equations~\eqref{eq:bal3}
%%
%\begin{equation}
%   \mathbf{J_{d \alpha}} = - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha =  - \left( D_\alpha \boldsymbol{\nabla} ( \rho Y_\alpha)
%   - \frac{D_\alpha}{\rho} \boldsymbol{\nabla} \rho \; ( \rho Y_\alpha) \right) \label{eq:expdfl}
%\end{equation}
%%
%which, substituted in the balance equation~\eqref{eq:bal3} leads to:
%%
%\begin{equation}
% \frac{\partial \rho Y_\alpha}{ \partial t} + \nabla \cdot \left(  \left[ \mathbf{v} + \frac{D_\alpha}{\rho} \boldsymbol{\nabla} \rho \right](\rho Y_\alpha)   \right) = \nabla \cdot \left( D_\alpha \boldsymbol{\nabla}  (\rho Y_\alpha ) \right) + \dot{m}_\alpha''' \; , \; \alpha=1,\dots,N \label{eq:bal4b}
%\end{equation}
%
%Next, we define the following implicit BE time discretization from time level $t$ (step $n$) to $t+\Delta t$ (step $n+1$) for the predictor step:
%%
%\begin{equation}
%\frac{ (\rho Y_\alpha)^{n+1} - (\rho Y_\alpha)^n}{\Delta t} = - \nabla \cdot \left(  \left[ \mathbf{v} + \frac{D_\alpha}{\rho} \boldsymbol{\nabla} \rho \right]^{n} (\rho Y_\alpha)^{n+1} \right) + \nabla \cdot \left( D_\alpha^n \boldsymbol{\nabla}  (\rho Y_\alpha )^{n+1} \right) + \dot{m}_\alpha'''^{n} \; , \; \alpha=1,\dots,N \label{eq:bal5}
%\end{equation}
%%
%where the advective and diffusive terms have been linearized, consistently with the time levels at which the velocity, diffusivity and mixture density are known, when solving the scalar transport-reaction equations within the predictor of the FDS scheme.
%The time discretization differs from a backward Euler discretization in the following:
%%
%\begin{itemize}
%  \item The velocity term $\mathbf{v} + \frac{D_\alpha}{\rho} \boldsymbol{\nabla} \rho$, is taken at step $n$ instead of step $n+1$.
%  \item Similarly, the diffusivity $D_\alpha$, is taken at step $n$.
%  \item The reactive term is known and taken at step $n$.
%\end{itemize}
%%
%
%A similar expression to~\eqref{eq:bal5} can be used when when implicitly solving on the corrector step with the trapezoidal rule, using the predicted (star) values to evaluate the factors on the different terms. That is, for example:


Consider implicit BE solution in cut-cell region for the first stage of SSPRK2 and trapezoidal rule applied in the second stage. Assume advective and diffusive terms have been linearized, consistently with the time levels at which the velocity, diffusivity and mixture density are known. The scheme defined in the previous section applied to species transport results in:
%
\begin{enumerate}
  \item Predictor:
%
\begin{equation}
\frac{ (\rho Y_\alpha)^{*} - (\rho Y_\alpha)^n}{\Delta t} = - \nabla \cdot \left(  \left[ \mathbf{v} +
\frac{D_\alpha}{\rho} \boldsymbol{\nabla} \rho \right]^{n} (\rho Y_\alpha)^{*} \right) +
\nabla \cdot \left( D_\alpha^n \boldsymbol{\nabla}  (\rho Y_\alpha )^{*} \right) + \dot{m}_\alpha'''^{n} \label{eq:pred}
\end{equation}
%
  \item Corrector:
%
\begin{equation}
\frac{ (\rho Y_\alpha)^{n+1} - (\rho Y_\alpha)^n}{\Delta t} = - \nabla \cdot \left(  \mathbf{v}' (\rho Y_\alpha) \right)^{n+1/2} + \nabla \cdot \left( D_\alpha \boldsymbol{\nabla}  (\rho Y_\alpha ) \right)^{n+1/2} + \dot{m}_\alpha'''^{n} \label{eq:corr}
\end{equation}
%
where the trapezoidal approximations to time level $t_{n+1/2}$ for advection and diffusion terms are:
%
\begin{eqnarray}
  \left(  \mathbf{v}'(\rho Y_\alpha) \right)^{n+1/2} & \simeq &
  \frac{1}{2} \left( \left[ \mathbf{v} + \frac{D_\alpha}{\rho} \boldsymbol{\nabla} \rho \right]^n (\rho Y_\alpha)^{n} +
                           \left[ \mathbf{v} + \frac{D_\alpha}{\rho} \boldsymbol{\nabla} \rho \right]^* (\rho Y_\alpha)^{n+1} \right)
                           \label{eq:ad12}  \\
  \left( D_\alpha \boldsymbol{\nabla}  (\rho Y_\alpha ) \right)^{n+1/2}   & \simeq &
  \frac{1}{2} \left(  D_\alpha^{n} \boldsymbol{\nabla}  (\rho Y_\alpha )^{n} +
                            D_\alpha^* \boldsymbol{\nabla}  (\rho Y_\alpha )^{n+1}  \right) \label{eq:df12}
\end{eqnarray}
%
\end{enumerate}
%
In the corrector discretization the lagged in time $\dot{m}_\alpha'''^{n}$ has been used, as defined in the new MOL (method of lines: advection diffusion discretized with SSPRK2, reaction uses forward Euler.) time integration scheme implemented in FDS. The mid-step values for the advection and diffusion factors~\eqref{eq:ad12}-\eqref{eq:df12} are approximated from the predicted unknowns, and values at step $n$. This stepping scheme combination is found to be second order time accurate in the numerical tests section.





\subsubsection{Cut-cell scheme for species transport implementation notes}


\subsection{Momentum time marching and Immersed Boundaries}

\subsubsection{Momentum-Pressure Coupling: IBM + CC}

As a first approximation, consider the Newtonian flow problem defined by the following set of partial differential equations:
%
\begin{eqnarray}
  \frac{\partial \mathbf{u}(\mathbf{x},t)}{\partial t} &=& - \left[ \mathbf{F}(\mathbf{u},\mathbf{x},t) + \boldsymbol{\nabla} H(\mathbf{x},t) \right] \; , \; \mathbf{x} \in \Omega - \sum{\Omega_i} \; , \; t \in \mathbb{R}_+ \label{eq:LowMachMom} \\
         \nabla \cdot \mathbf{u} (\mathbf{x},t) & = & \left(\nabla \cdot \mathbf{u} \right)^{th} \label{eq:LowMachDiv}
\end{eqnarray}
%
where equation~\eqref{eq:LowMachMom} in the momentum equation, $\mathbf{u}(\mathbf{x},t)$ is the spatial velocity field, $\mathbf{F}(\mathbf{u},\mathbf{x},t)$ is a vector containing convective, diffusive and possibly other force terms, and $H(\mathbf{x},t)$ is a potential scalar field (physically the head field in this case, commonly called pressure). For sake of argument here, to represent the low Mach approximation employed in FDS, it suffices to consider a specified divergence field $\left(\nabla \cdot \mathbf{u} \right)^{th} (\mathbf{x},t)$ (thermodynamic divergence). This divergence field in the thermally buoyant flow model used in FDS is a proxy for the energy equation.
The domain $\Omega - \sum{\Omega_i}$ represents the fluid region, and boundary conditions are prescribed for $\mathbf{u}(\mathbf{x},t)$ on $\partial \Omega,\partial \Omega_1,...,\partial \Omega_{nbods}$.

Classical fractional step methods for time integration of incompressible or low Mach flow are based on two operations: First, momentum transport to obtain intermediate velocities, and second, projection of velocities into divergence free space. Consider the Forward Euler update of the governing equations from $t_n$ to $t_{n+1}=t_n + \Delta t$ of the form: Given $ \mathbf{u}^n=\mathbf{u}(\mathbf{x},t_n)$, $\nabla \cdot \mathbf{u}^{n+1} = \left( \nabla \cdot \mathbf{u}^{n+1} \right)^{th}$ known
%
\begin{eqnarray}
  \frac{\mathbf{u}^{n+1}-\mathbf{u}^{n}}{\Delta t} &=& - \left[ \mathbf{F}^n +  \boldsymbol{\nabla} H^n \right] \label{eq:LowMachMomEu}\\
  \nabla \cdot \mathbf{u}^{n+1} &=& \left( \nabla \cdot \mathbf{u}^{n+1} \right)^{th} \label{eq:LowMachDivEu}
\end{eqnarray}
%
where as time has been discretized, $\mathbf{u}^{n+1}$ represents a numerical approximation to the solution in~\eqref{eq:LowMachMom}-\eqref{eq:LowMachDiv} at time $t_{n+1}$. As the potential field $H(\mathbf{x},t)$ does not have a time evolution equation, it is assumed responsible of enforcing the divergence condition and used on the projection step. Taking the divergence of equation~\eqref{eq:LowMachMomEu} and considering  the constraint~\eqref{eq:LowMachDivEu}, the two steps of the method are
%
\begin{enumerate}
  \item Solve Poisson equation for $H^n$:

\begin{equation}
   \nabla \cdot \boldsymbol{\nabla} H^n = - \left[ \frac{\left( \nabla \cdot \mathbf{u}^{n+1} \right)^{th} - \nabla \cdot \mathbf{u}^{n}}{\Delta t} \right] - \nabla \cdot \mathbf{F}^n \label{it:FSPoisson}
\end{equation}

  \item Obtain final velocity for step:

  \begin{equation}
     \mathbf{u}^{n+1} = \mathbf{u}^{n} - \Delta t \left[ \mathbf{F}^n +  \boldsymbol{\nabla} H^n \right] \label{it:FSProject}
   \end{equation}

   The term $\hat{\mathbf{u}}^{n+1}=\mathbf{u}^{n} - \Delta t \mathbf{F}^n$ is known as intermediate velocity, and is a non matching divergence approximation to $\mathbf{u}^{n+1}$ (i.e. $\nabla \cdot \hat{\mathbf{u}}^{n+1} \neq \left( \nabla \cdot \mathbf{u}^{n+1} \right)^{th}$).
\end{enumerate}
%
Although in the original problem~\eqref{eq:LowMachMom}-\eqref{eq:LowMachDiv} no boundary condition is required for $H(\mathbf{x},t)$, a consequence of the projection scheme is that boundary conditions are required on the Poisson equation of step~\eqref{it:FSPoisson}. For explicit methods and stationary \textit{solid} boundaries, the corresponding boundary condition is \textit{homogeneous} Neumann for $H^n$, $\partial H^n / \partial x_n =0$ in $\partial \Omega,\partial \Omega_1,...,\partial \Omega_{nbods}$ (i.e.~\cite{Perot:1993}).

%Indeed, when the velocity boundary condition is applied to the intermediate velocity $u^*=\mathbf{u}^{n} - \Delta t\mathbf{F}^n$, mass conservation on step~\ref{it:FSProject} requires that at the boundary $(\mathbf{u}^{n+1} -u^*) \cdot \mathbf{\hat{n}} =0$, from which the previously defined condition for $H^n$ arises.

%Another approach to projection schemes is based on the method of lines~\cite{schiesser1991numerical}. It assumes initially a spatial discretization of equations~\eqref{eq:insMom}-\eqref{eq:insDiv}, and as second step temporal discretization of the resulting system of ODEs, employing discrete projection.
%
%We consider first a rectangular domain with no immersed boundaries. When using a regular uniform staggered mesh and a central finite difference method with second order spatial accuracy for discretization the spatially discrete version of~\eqref{eq:insMom}-\eqref{eq:insDiv} is:
%%
%\begin{eqnarray}
%  [\mathbf{M}_{FD}] \{\mathbf{\dot{U}}\} &=& - \left( [\mathbf{C}_{FD(\{\mathbf{U}\})}] \{\mathbf{U}\} +
%  [\mathbf{L}_{FD}] \{\mathbf{U}\} + [\mathbf{G}_{FD}] \{\mathbf{H}\} \right) + \{\mathbf{f}_{bcFD}\} \label{eq:momspdisc}\\
%                                                                         \left[\mathbf{D}_{FD}\right]     \{\mathbf{U}\} & = & \{\mathbf{0}\} + \{\mathbf{d}_{bcFD}\}      \label{eq:divspdisc}
%\end{eqnarray}
%%
%here $[\mathbf{M}_{FD}]=[\mathbf{I}]$  of size $n_u \times n_u$ where $n_u$ is the number of velocity unknowns, is the mass matrix for the finite difference discretization, and $[\mathbf{C}_{FD(\{\mathbf{U}\})}],[\mathbf{L}_{FD}]$ are momentum convection and diffusion matrices. The scalar gradient matrix $[\mathbf{G}_{FD}]$ is of size $n_u \times n_h$ where $n_h$ is the number of scalar $H$ unknows, and the divergence matrix $\left[\mathbf{D}_{FD}\right]$ is of size   $n_h \times n_u$. At this point boundary conditions have been applied modifying the named matrices and adding vectors $\{\mathbf{f}_{bcFD}\},\{\mathbf{d}_{bcFD}\}$ in case of nonhomogeneous conditions.
%Using a Forward Euler time discretization for equations~\eqref{eq:momspdisc}-\eqref{eq:divspdisc}, representative of explicit schemes we obtain
%%
%\begin{equation}
%  \left[ \begin{array}{cc} \mathbf{A}_{FD} & \mathbf{G}_{FD} \\
%                                     \mathbf{D}_{FD}& \mathbf{0} \end{array} \right]
%  \left\{ \begin{array}{c}  \mathbf{U}^{n+1}  \\  \mathbf{H}^n \end{array} \right\}   =
%  \left\{ \begin{array}{c}  \mathbf{A}_{FD}\mathbf{U}^{n} - \mathbf{F}^{n}_{FD}  \\ \mathbf{0}  \end{array} \right\} +
%  \left\{ \begin{array}{c} \mathbf{f}_{bcFD}  \\ \mathbf{d}_{bcFD}   \end{array} \right\} \label{eq:saddle}
%\end{equation}
%%
%where $[\mathbf{A}_{FD}] =1/\Delta t [\mathbf{M}_{FD}]$ is invertible, $\{\mathbf{F}^{n}_{FD}\}=[\mathbf{C}_{FD(\{\mathbf{U^n}\})}] \{\mathbf{U^n}\} +[\mathbf{L}_{FD}] \{\mathbf{U^n}\}$. The matrix on the right hand side is block LU decomposed as~\cite{Perot199351}
%%
%\begin{equation}
%   \left[ \begin{array}{cc} \mathbf{A}_{FD} & \mathbf{G}_{FD} \\
%                                     \mathbf{D}_{FD}& \mathbf{0} \end{array} \right]  \simeq
%  \left[ \begin{array}{cc} \mathbf{A}_{FD} & \mathbf{0} \\
%                                     \mathbf{D}_{FD}& - \mathbf{D}_{FD} \mathbf{B}_{FD} \mathbf{G}_{FD} \end{array} \right]
%  \left[ \begin{array}{cc} \mathbf{I} &  \mathbf{B}_{FD} \mathbf{G}_{FD} \\
%                                     \mathbf{0}& \mathbf{I} \end{array} \right]    \label{eq:fslu}
%\end{equation}
%%
%For the previous factorization to be exact $[\mathbf{B}_{FD}]=[\mathbf{A}_{FD}]^{-1}$ (i.e. Uzawa method~\cite{Chang2002183}). Classic fractional step schemes use a first order approximation for  $[\mathbf{B}]=\Delta t [\mathbf{I}]$, which in general results in loss of temporal accuracy for velocities when high order implicit methods are used for time integration~\cite{Perot199351,Chang2002183}. For explicit time integration, and finite difference or difference finite volume spatial discretizations, matrix $[\mathbf{M}]$ is lumped (i.e. \textit{diagonal}) as is the case of lumped element finite volume discretizations. Then,  $[\mathbf{B}]=[\mathbf{A}]^{-1}=\Delta t [\mathbf{M}]^{-1}$ is easily obtained. Indeed, $[\mathbf{B}_{FD}]=\Delta t [\mathbf{I}]$ and $[\mathbf{B}_{FV}]=\Delta t / V_{cell} [\mathbf{I}]$ for finite differences and uniform grid difference finite volume respectively ($V_{cell}$ uniform cell volume).
%
%In terms of intermediate velocities $\{ \mathbf{U}^*\}$ the LU decomposed problem~\eqref{eq:saddle}-\eqref{eq:fslu} is divided in
%%
%\begin{eqnarray}
%  [\mathbf{A}_{FD}] \{ \mathbf{U}^*\} &=& [\mathbf{A}_{FD}] \{\mathbf{U}^{n}\} - \{\mathbf{F}^{n}_{FD}\} + \{\mathbf{f}_{bcFD}\} \\
%  \left[ \mathbf{D}_{FD} \mathbf{B}_{FD} \mathbf{G}_{FD} \right] \{\mathbf{H}^n\} &=&  [\mathbf{D}_{FD}] \{ \mathbf{U}^*\} - \{\mathbf{d}_{bcFD}\} \\
%  \{\mathbf{U}^{n+1}\} &=& \{ \mathbf{U}^*\} - [\mathbf{B}_{FD}] [\mathbf{G}_{FD}] \{\mathbf{H}^n\}
%\end{eqnarray}
%Considering homogenous boundary conditions, the zero discrete divergence of $\{\mathbf{U}^{n}\}$ and $[\mathbf{M}_{FD}]=[\mathbf{I}]$ we recover the FD spatially discretized version of steps~\ref{it:FSPoisson}-\ref{it:FSProject}:
%%
%\begin{eqnarray}
%  \left[ \mathbf{D}_{FD} \mathbf{G}_{FD} \right] \{\mathbf{H}^n\} &=& - [\mathbf{D}_{FD}] \{\mathbf{F}^{n}_{FD}\} \\
%  \{\mathbf{U}^{n+1}\} &=& \{\mathbf{U}^{n}\} - \Delta t \left( \{\mathbf{F}^{n}_{FD}\} + [\mathbf{G}_{FD}] \{\mathbf{H}^n\} \right)
%\end{eqnarray}
%%
The next component of our scheme is used to approximate the no slip boundary condition for immersed solid boundaries. We employ a direct forcing immersed boundary method for the momentum equations~\cite{Fadlun:2000}. In order to do this, a force field is devised on the discrete momentum equations on grid faces crossed by the immersed surfaces, to approximate the no-slip boundary condition on said surfaces. We describe some preliminary concepts before introducing our IBM+CC momentum-pressure coupling scheme.

\subsubsection{Simple Planar Velocity interpolation:}

Immersed boundary forcing requires estimates for velocities at time step $n+1$ on \texttt{INBOUNDARY} and \texttt{GASPHASE} cut-faces. For the former, the face velocities are assumed known at step $n+1$, given by the known (or predicted) motion of the immersed bodies (i.e. horizontal velocity component $u_{B}^{n+1}$ in immersed boundary point $B$ defined by $\mathbf{x}_B$ in figure~\ref{Fig:IBCart}). For collocated velocities on  \texttt{GASPHASE} cut-face centroids interpolation is required. For linear velocity interpolation from a boundary point $B$ and a $\{1,...,ne\}$ set of \textit{fluid} points to forcing point $d$, the target $x$ velocity component $u_d$ on gas phase cut-face $\mathbf{3}$ of figure~\ref{Fig:IBCart} is defined by the following interpolation formula:
%
\begin{equation}
   u_d = \sum_{e=0}^{ne} \phi (\mathbf{x}_d - \mathbf{x}_e) u_e \label{eq:interpu}
\end{equation}
%
where $u_0=u_B$ and $\mathbf{x}_0=\mathbf{x}_B$, and the rest of points in the interpolation stencil are \textit{fluid} points.
Note: Fluid points are in this case, velocity points in regular fluid faces which are not being forced by the immersed boundary scheme.

We describe a simple interpolation method for $u_d$ in the two dimensional case of figure~\ref{Fig:IBCart}, using a variant of the normal probe approach~\cite{Balaras:2004}. In this case, boundary point $B$ is defined as the closest point to forcing point $d$, and is either a point within a triangle whose normal $x_n$ direction is parallel to segment $\overline{B\mbox{-} d}$, or an edge or vertex point defining a direction $x_n$ into the fluid domain.
Define a point $\overline{e_x}$ by the intersection of line $\mathbf{l}(x_n)=\mathbf{x}_B+ x_n \mathbf{\hat{n}}$ and the first grid line containing two fluid points $e=1,2$ which surround said intersection point. The grid line in 2D chosen is the one defining a Cartesian plane with normal most aligned with $\mathbf{\hat{n}}$. See figure~\ref{Fig:IBCart}.
The procedure to obtain an approximation for $\hat{u}_d$ at step $n+1$ is as follows:

%
\begin{enumerate}

  \item Approximate velocity at $t_{n+1}$ on \textit{fluid points} $e=1,2$ as: $\hat{u}_e^{n+1}= u_e^{n} -
  \Delta t \left[ F_e^n + (\partial H^{n,k-1} / \partial x)_e \right]$. Here, the head gradients $(\partial H^{n,k-1} / \partial x)_e$ are assumed known.
  Use linear interpolation on grid line $L$ to estimate $\hat{u}_{ex}^{n+1}$:
  \begin{equation}
     \hat{u}_{ex}^{n+1} = C_1 \hat{u}_1^{n+1} + C_2 \hat{u}_2^{n+1} \label{eq:intex}
  \end{equation}
  where $C_1=||\mathbf{x}_2-\mathbf{x}_{ex}||/||\mathbf{x}_2-\mathbf{x}_{1}||$, $C_2=||\mathbf{x}_1-\mathbf{x}_{ex}||/||\mathbf{x}_2-\mathbf{x}_{1}||$. Note that $C_1+C_2=1$.

  \item Surface velocity $u_B^{n+1}$ is known.

  \item Apply linear interpolation for cut-face velocity approximation at $n+1$:
  \begin{equation}
     \hat{u}_d^{n+1}=D_B u_B^{n+1} + D_{ex} \hat{u}_{ex}^{n+1} \label{eq:intd}
  \end{equation}
   where $D_B=||\mathbf{x}_{ex}-\mathbf{x}_{d}||/||\mathbf{x}_{ex}-\mathbf{x}_{B}||$, $D_{ex}=||\mathbf{x}_{B}-\mathbf{x}_{d}||/||\mathbf{x}_{ex}-\mathbf{x}_{B}||$.

  Note that $D_B+D_{ex}=1$.  Using equation~\eqref{eq:intex} in equation~\eqref{eq:intd} we have:
  \begin{equation}
      \hat{u}_d^{n+1}=D_B \; u_B^{n+1} + D_{ex} C_1 \; \hat{u}_1^{n+1} + D_{ex} C_2 \; \hat{u}_2^{n+1}  \label{eq:interpu2}
  \end{equation}
  where by the previous the sum of the interpolation coefficients $D_{B} + D_{ex} C_1 + D_{ex} C_2 = 1$. Equation~\eqref{eq:interpu2} is a specialized version of the general expression~\eqref{eq:interpu} for the interpolation scheme described.

\end{enumerate}
%
Similar expressions are used for velocity components on other directions (i.e. $\hat{v}_g^{n+1}$ in figure~\ref{Fig:IBCart}). This interpolation procedure can be naturally be expanded to three dimensions, where a set of velocity fluid points in an outer Cartesian grid plane can be used, such that a two dimensional interpolation on the plane is performed.

% figure
\begin{figure}[h]
      \centering
      \includegraphics[trim = 30mm 22mm 40mm 22mm, clip,
       width=0.75\linewidth]{../../../fig/fds/Geometry_Figures/CutCellsIBSketchNormal.png}
      \caption{Velocity interpolation along normal through cut-face centroid $d$: \texttt{GASPHASE} cut-cells,  corresponding \texttt{GASPHASE} regular and cut-faces (in red), and underlying Eulerian staggered mesh are shown. Velocity points used for linear IB interpolation to cut-face centroids (in red) are shown. }
	\label{Fig:IBCart}
\end{figure}


\subsubsection{Dynamic Velocity estimation for wall modeled IBM:}

In this section we describe a scheme to dynamically define the value of $\hat{u}_d^{n+1}$ making use of an equilibrium wall model solution and simplified boundary layer equations along a local streamwise direction.

\bigskip
\noindent
\textbf{Boundary conforming grids:}
\bigskip

In order to define this scheme, we first look at the case of wall modeling with a simple wall function (\textit{log law}) for a boundary conforming mesh. Consider the case depicted in figure X, where a velocity component $\hat{u}_e^n$ is defined in a face adjacent to the solid boundary. The vertex point on the wall for this face is point $B$. A simple model representation in this case is given by a function $u^+ = f(y^+)$:
%
\begin{eqnarray}
   u^+ = y^+ \; \; , \; \; y^+ < 12 \\
   u^+ = \frac{1}{\kappa} ln(y^+) + Br \; \; , \; \; y^+ > 12
\end{eqnarray}
%
and the relationships:
%
\begin{equation}
   u_\tau = \sqrt{\frac{\tau_w}{\rho}}; \; \; y^+= \frac{u_\tau y}{\nu}; \; \; u^+=\frac{u}{u_\tau}
\end{equation}
%
which approximate the mean solution of an equilibrium boundary layer problem. Here $\tau_w$ is the stress at the wall point $B$; $\rho, \nu$ are the density and molecular kinematic viscosity, and $\kappa \simeq 0.41, Br \simeq 5.2$ are the log law constants.

Through the expression for the wall stress $\tau_w = \rho \nu \left( \frac{d u}{d y} \right)_{y=0}$ and a known value $\hat{u}_e^n$, a fixed iteration can be applied to compute the corresponding $y^{+}_e$. Start by:
%
\begin{itemize}
   \item compute $\tau_w^0 = \rho \nu \frac{2 |\hat{u}_e^n|}{\Delta y}$
   \item compute $u_\tau^0 = \sqrt{\tau_w^0 / \rho}$
   \item compute $\delta_\nu^0 = \frac{\nu}{u_\tau^0}$
   \item compute $y^{+,0} = \frac{\Delta y}{2 \delta_\nu^0}$
\end{itemize}
%
If a direct numerical simulation (DNS) is performed, this is all we need to compute the wall stress (its definition). If we perform a wall modeled large eddy simulation (LES), for $k=1,2,..$ up to desired convergence do:
%
\begin{enumerate}
   \item compute $y^{+,k} = \frac{\Delta y}{2 \delta_\nu^{k-1}}$
   \item if $y^+ < 12$ : compute $\tau_w^k = \rho (\hat{u}_e^n / y^{+,k})^2$; $u_\tau^k = \sqrt{\tau_w^k/\rho}$ %; $\left( \frac{\partial u}{\partial y} \right)_w^k = \frac{2 |\hat{u}_e^n|}{\Delta y}$
   \item else : compute   $\tau_w^k = \rho \left( \frac{\hat{u}_e^n}{\ln (y^{+,k})/\kappa+Br} \right)^2$; $u_\tau^k = \sqrt{\tau_w^k/\rho}$ %; $\left( \frac{\partial u}{\partial y} \right)_w^k = \frac{2 u_\tau}{\kappa \Delta y}$
   \item compute $\delta_\nu^{k} = \frac{\nu}{u_\tau^k}$
\end{enumerate}
%
%The final value of $\left( \frac{\partial U}{\partial y} \right)_w^k$ can be used to compute a slip factor $S_F$, that would define the ghost cell value of velocity, in order to compute the right stress from the centered finite difference definition. This is done in the following form:
%%
%\begin{equation}
%   \left( \frac{\partial u}{\partial y} \right)_w = \frac{\hat{u}_e^n - \hat{u}_g}{\Delta y} = (1-S_F) \frac{\hat{u}_e^n}{\Delta y}
%\end{equation}
%%
%from which the slip factor can be computed as $S_F = 1 -  \left( \frac{\partial u}{\partial y} \right)_w \frac{\Delta y}{|\hat{u}_e^n|}$.
%
Once the stress at the wall $\tau_w$ is computed, it is used in the right hand side definition of the momentum equations to evolve the velocity flow variable $\hat{u}_e^n \rightarrow \hat{u}_e^{n+1}$. Note the velocity on this face point off the wall is integrated in time in the same manner as the rest of the discrete velocities, having to obey CFL and viscous constraints for stability.

% figure: Boundary conforming setup:



\bigskip
\noindent
\textbf{Momentum equations along a local streamwise direction}
\bigskip

The simplified equations used in a streamwise - wall normal local coordinate system are:
%
\begin{eqnarray}
   \frac{\partial u_{\hat{\mathbf{s}}}}{\partial s} +  \frac{\partial u_{\hat{\mathbf{n}}}}{\partial n} &=& \left( \nabla \cdot \mathbf{u} \right)^{th}  \label{eq:blcont} \\
    \frac{\partial u_{\hat{\mathbf{s}}}}{\partial t} + u_{\hat{\mathbf{s}}} \frac{\partial u_{\hat{\mathbf{s}}}}{\partial s} +
    u_{\hat{\mathbf{n}}} \frac{\partial u_{\hat{\mathbf{s}}}}{\partial n} &=& - \frac{1}{\rho} \left( \frac{\partial p}{\partial s} +
    \frac{\partial \tau_{sn}}{\partial n} - (\rho -\rho_0) g_s \right) \label{eq:blmom}
\end{eqnarray}
%
where $\hat{\mathbf{s}},\hat{\mathbf{n}}$ are the streamwise and normal unit vectors, and $u_{\hat{\mathbf{s}}},u_{\hat{\mathbf{n}}}$ are the corresponding velocity components. Also, $p$ is the pressure field, $\tau_{sn}$ is the shear stress in local coordinates and $g_s$ is the gravity acceleration along the streamwise direction.

\bigskip
\noindent
\textbf{Geometry not conforming to fluid mesh:}
\bigskip

For cases where the geometry does not conform to the fluid mesh, as depicted in figure X1, the sequence of operations to define an immersed boundary forcing term in gas phase cut faces is the following:
%
\begin{itemize}
   \item Define normal probe components for gas cut-face centroid $d$, in particular, boundary point B and external point $e_x$.
   \item Interpolate known velocities, velocity gradients and cell centered variables to point $e_x$.
   \item Define a local coordinate system represented by the normal $\hat{\mathbf{n}}$, streamwise $\hat{\mathbf{s}}$ and orthonormal $\hat{\mathbf{t}}$ unit vectors.
   \item Apply wall model described before to estimate the known streamwise velocity $u_{\hat{\mathbf{s}},d}^n$ value at point $d$ and wall stress $\tau_w$ in this local coordinate system.
   \item Estimate pressure gradient and gravity force along streamwise direction.
   \item Using streamwise-normal boundary layer equations: estimate velocity in the normal direction $u_{\hat{\mathbf{n}},d}^{n+1}$, and perform a one step time integration of the streamwise boundary layer equations to obtain $u_{\hat{\mathbf{s}},d}^{n+1}$.
   \item Recompute the component $\hat{u}_d^{n+1}$ on the gas cut-face.
\end{itemize}
%
We note the previous steps can be gathered in four general tasks: interpolation, wall function estimation of wall stress, integration in local coordinate system, and recomputation of $\hat{u}_d^{n+1}$ on the gas cut-face (inertial reference frame).

\bigskip
\noindent
\textbf{Interpolation of variables to external point:}
\bigskip

The first task, interpolation, requires the following operations \textit{for each} gas cut-face:
%
\begin{itemize}
   \item Define external point $e_x$ in the normal direction crossing through the gas cut-face centroid $d$.
   \item Use an interpolation stencil composed of fluid points to interpolate velocity components and gradients from staggered locations to $e_x$ (i.e. equation~\eqref{eq:interpu}).
   \item Use an interpolation stencil to interpolate total viscosity to point $e_x$ from Cartesian cell centroids. Similar operation to obtain a pressure gradient.
\end{itemize}
%

We describe the trilinear interpolation of one velocity component to the external point $e_x$, masked to a stencil of admissible fluid points.
The first step is to locate the point $e_x$ within the cuboid space defined by  eight $u$ velocity component collocation points $\mathbf{x}_k$, $k=1,..,8$, in the fluid grid. The interpolation coefficients from these to $e_x$ are noted as $\phi_k(\mathbf{x}_{ex}-\mathbf{x}_k)$. Consider the $\mathbf{x}_k$ collocation points defining the vertices of a box with coordinates $[x_{lo},x_{hi}] \times [ y_{lo},y_{hi}] \times [z_{lo} , z_{hi}]$, and the normalized coordinates $\xi,\eta,\zeta$:
%
\begin{equation}
   \xi = \frac{x_{ex}-x_{lo}}{x_{hi}-x_{lo}}, \; \; \;  \eta = \frac{y_{ex}-y_{lo}}{y_{hi}-y_{lo}}, \; \; \; \zeta = \frac{z_{ex}-z_{lo}}{z_{hi}-z_{lo}}
\end{equation}
%
Then, $\phi_k(\mathbf{x}_{ex}-\mathbf{x}_k)$ takes the form:
%
\begin{eqnarray}
    \phi_k(\mathbf{x}_{ex}-\mathbf{x}_k) =
    \left[ \frac{x_{k}-x_{lo}}{x_{hi}-x_{lo}} \xi + \frac{x_{hi}-x_{k}}{x_{hi}-x_{lo}} (1-\xi) \right] \times \nonumber \\
    \left[ \frac{y_{k}-y_{lo}}{y_{hi}-y_{lo}} \eta + \frac{y_{hi}-y_{k}}{y_{hi}-y_{lo}} (1-\eta) \right] \times \nonumber \\
    \left[ \frac{z_{k}-z_{lo}}{z_{hi}-z_{lo}} \zeta + \frac{z_{hi}-z_{k}}{z_{hi}-z_{lo}} (1-\zeta) \right] \; \; \; \; \;
\end{eqnarray}
%

Depending on the distance $\delta_{ex}$ from the boundary to point $e_x$, some of these $k$ points might lay in Cartesian faces containing gas cut-faces, which will be eventually also forced. To avoid this dependence, a mask is defined in the interpolation, such that the final interpolation coefficients used are non-zero only in fluid points. We require the nonzero masked interpolation coefficients to be positive and add to unity. Ideally we would want them to preserve the accuracy of the original interpolation.
We rescale the interpolation coefficients on fluid points, dividing them by their sum. Also, we set to zero (mask) the interpolation coefficients related to inadmissible points.

Other unstructured interpolation schemes can also be implemented defined only on fluid points (i.e. weighted least squares, inverse distance method).


\bigskip
\noindent
\textbf{Local streamwise - wall normal coordinate system velocity components:}
\bigskip

Having the wall normal direction $\hat{\mathbf{n}}$, an orthonormal direction can be computed from:
%
\begin{equation}
    \hat{\mathbf{t}} = \frac{\hat{\mathbf{n}} \times \mathbf{u}_{ex}^n}{\| \hat{\mathbf{n}} \times \mathbf{u}_{ex}^n \| }
\end{equation}
%
where $\mathbf{u}_{ex}^n=\mathbf{u}_{ex,abs}^n-\mathbf{u}_{B}^n$ is the known relative velocity of the fluid respect to the boundary, and the absolute velocity $\mathbf{u}_{ex,abs}^n$ was interpolated as described previously.
The streamwise velocity component is defined by the product $\hat{\mathbf{s}} = \hat{\mathbf{t}} \times \hat{\mathbf{n}}$ and streamwise and normal velocity components  $u_{\hat{\mathbf{s}},ex}^n= \mathbf{u}_{ex}^n \cdot \hat{\mathbf{s}}$, $u_{\hat{\mathbf{n}},ex}^n= \mathbf{u}_{ex}^n \cdot \hat{\mathbf{n}}$. The local component in the $\hat{\mathbf{t}}$ direction is zero by construction.

\bigskip
\noindent
\textbf{Wall function estimation of wall stress $\tau_w$:}

With the local streamwise velocity component, as well as fluid properties $\rho, \mu$ and the distance $\delta_{ex}$ from boundary to the $ex$ point, we are equipped to compute the wall function iteration described for boundary conforming grids. From this iteration, values of $\tau_w$, $y^+$ and $u_\tau$ and $u_{\hat{\mathbf{s}},d}^n$ can be computed. The last variable is the velocity at time step $n$ in point $d$, and along the streamwise direction.

\bigskip
\noindent
\textbf{Cut-face centroid velocity $\hat{u}_d^{n+1}$ update from local momentum equations}
\bigskip

We discretize equations \eqref{eq:blcont}-\eqref{eq:blmom} in the following form:
%
\begin{eqnarray}
   u_{\hat{\mathbf{n}},d}^{n+1} &=& \mathbf{u}_{B}^n \cdot \hat{\mathbf{n}} + \delta_d \left(  \left( \nabla \cdot \mathbf{u} \right)^{th}
   - \frac{\partial u_{\hat{\mathbf{s}}} }{\partial s} |_{d}^n \right) \label{eq:un} \\
   RHS_{u_{\hat{\mathbf{s}},d}} &=&  u_{\hat{\mathbf{s}},d}^n  \frac{\partial u_{\hat{\mathbf{s}}} }{\partial s} |_{d}^n  +
   u_{\hat{\mathbf{n}},d}^{n+1} \frac{\partial u_{\hat{\mathbf{s}}} }{\partial n} |_{d}^n +
   \frac{1}{\rho} \left[  \frac{\partial p}{\partial s} |_d^n + \frac{\tau_{sn} |_{ex}^n - \tau_w}{\delta_{ex}} - (\rho-\rho_0) g_s \right] \label{eq:rhsus} \\
   u_{\hat{\mathbf{s}},d}^{n+1} &=& u_{\hat{\mathbf{s}},d}^{n}+\mathbf{u}_{B}^n \cdot \hat{\mathbf{s}} - \Delta t RHS_{u_{\hat{\mathbf{s}},d}} \label{eq:us}
\end{eqnarray}
%
In equations~\eqref{eq:un}-\eqref{eq:us} we compute:
%
\begin{eqnarray}
   \frac{\partial u_{\hat{\mathbf{s}}} }{\partial s} |_{d}^n &\simeq& \frac{ \delta_d}{ \delta_{ex}} \frac{\partial u_{\hat{\mathbf{s}}} }{\partial s} |_{ex}^n \\
    \frac{\partial u_{\hat{\mathbf{s}}} }{\partial n} |_{d}^n &\simeq& \frac{ \delta_d}{ \delta_{ex}} \frac{\partial u_{\hat{\mathbf{s}}} }{\partial n} |_{ex}^n
\end{eqnarray}
%
that is, a linear variation of the streamline velocity component $s$ derivative along the wall normal direction has been assumed.
The term $\frac{\partial u_{\hat{\mathbf{s}}} }{\partial s} |_{ex}^n$ on the previous is computed from the expression:
%
\begin{eqnarray}
   \frac{\partial u_{\hat{\mathbf{s}}} }{\partial s} |_{ex}^n &=& \hat{\mathbf{s}}^T \left[ \nabla \mathbf{u}_{ex}^n \right]  \hat{\mathbf{s}} \\
   \left[ \nabla \mathbf{u}_{ex}^n \right] &=& \begin{bmatrix}
   \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} \\
   \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\
   \frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z} \\
   \end{bmatrix}_{ex}^n
\end{eqnarray}
%
where the derivatives of velocity components at the external $e_x$ point are computed from derivatives of the interpolation functions. Also,  $\frac{\partial u_{\hat{\mathbf{s}}} }{\partial n} |_{ex}^n=\hat{\mathbf{n}}^T \left[ \nabla \mathbf{u}_{ex}^n \right]  \hat{\mathbf{s}}$.

Further, $\frac{\partial p}{\partial s}  |_d^n \simeq \frac{\partial p}{\partial s} |_{ex}^n = (\nabla p)_{ex}^n \cdot \hat{\mathbf{s}}$,
and, being $[ \tau ]_{ex}^n$ the stress tensor at $e_x$ evaluated from interpolated velocity derivatives and total viscosity, the shear component of the rotated stress tensor $\tau_{sn} |_{ex}^n = \hat{\mathbf{s}}^T [ \tau ]_{ex}^n \hat{\mathbf{n}}$. Finally, on the buoyancy term, $\rho$, $\rho_0$ are interpolated to the $d$ centroid position,  and $g_s = \mathbf{g} \cdot \hat{\mathbf{s}}$.


The target velocity on the cut-face centroid, assumed in the $i$ Cartesian direction is:
%
\begin{equation}
   u_d^{n+1} = u_{\hat{\mathbf{n}},d}^{n+1} (\hat{\mathbf{n}} \cdot \hat{\mathbf{e}}_i) + u_{\hat{\mathbf{s}},d}^{n+1} (\hat{\mathbf{s}} \cdot \hat{\mathbf{e}}_i) + (\mathbf{u}_{B}^{n+1} \cdot \hat{\mathbf{t}})(\hat{\mathbf{t}} \cdot \hat{\mathbf{e}}_i)
\end{equation}
%
where $\hat{\mathbf{e}}_i$ is the Cartesian versor in direction $i$.

\subsubsection{Equivalence of discrete divergence integrals among cut-cells and underlying Eulerian cells:}

Consider Eulerian cell $i,j$ in figure~\ref{Fig:IBCart}, containing a \texttt{GASPHASE} cut-cell $ii$ with centroid $\mathbf{x}_{ii}$. The integral of the Eulerian cell velocity divergence can be divided as:
%
\begin{equation}
   \int_{\Omega_{i,j}}{ \nabla \cdot \mathbf{u}} \; d\Omega_{i,j} = \int_{\Omega_{ii}}{ \nabla \cdot \mathbf{u}} \; d\Omega_{ii} +
                                                                                                     \int_{\Omega_{kk}}{ \nabla \cdot \mathbf{u}} \; d\Omega_{kk}
\end{equation}
%
Where $kk$ refers to the solid region of cell $(i,j)$. Assume the solid region behaves as a rigid solid material, and therefore $\int_{\Omega_{kk}}{ \nabla \cdot \mathbf{u}} \; d\Omega_{kk} = 0$.
Referring to figure~\ref{Fig:IBCart}, and using the divergence theorem for each of the previous integrals, we have for cell $i,j$
%
\begin{eqnarray}
  \int_{\Omega_{i,j}}{ \nabla \cdot \mathbf{u}} \; d\Omega_{i,j} &=& \sum_{k \in \partial \Omega_{i,j}}{(\mathbf{u} \cdot \mathbf{\hat{n}})_k A_k} = (u_E-u_W) \Delta y + (v_N-v_S) \Delta x \label{eq:intij}\\
  \int_{\Omega_{ii}}{ \nabla \cdot \mathbf{u}} \; d\Omega_{ii}  &=& u_d \alpha_y \Delta y + v_d \alpha_x \Delta x -
  \gamma_{B1} A_1 - \gamma_{B2} A_2   \label{eq:intii}
\end{eqnarray}
%
where $\alpha_y=A_3/\Delta y$, $\alpha_x=A_4/\Delta x$ are the corresponding  \texttt{GASPHASE} area factors. The components $ \gamma_{B1}, \gamma_{B2}$ are the average normal velocities into the gas region for  \texttt{INBOUNDARY} cut-faces $\mathbf{1},\mathbf{2}$. Noting that for a stationary body all face velocities on solid region are zero (i.e. $u_W=v_S=\gamma_{B1}=\gamma_{B2}=0$) and equating~\eqref{eq:intii} with~\eqref{eq:intij} we find that
%
\begin{eqnarray}
  u_E = \alpha_y \: u_d  \label{eq:uE} \\
  v_N  = \alpha_x \: v_d  \label{eq:vN}
\end{eqnarray}
%
which express the velocities in the Eulerian cell faces in terms of velocities on fluid and solid cut-faces for matching flux and discrete divergence integrals. Therefore, knowing the velocities on the fluid (and solid) portion of cut-faces we can define the flux matched corresponding velocities on the underlying Eulerian grid face, such that the discrete divergence integral on this zone is preserved.
These expressions allow for the transfer of velocity among cut-faces and underlying Eulerian grid faces in a discrete divergence preserving manner.
Knowing the geometry parameters for all cut-faces, this equivalence can be extended for all cut-cells in a straightforward manner.

\subsubsection{IBM + CC coupling using discrete divergence preserving velocity transfer and underlying structured mesh}

This method uses a combination of finite difference discretization on regular \texttt{GASPHASE} and cut-cell underlying Eulerian cells, where both face velocities and $H$ at cell centers are defined. Transfer of velocities from \texttt{GASPHASE} cut-face centroids to Eulerian faces is done on a flux conserving manner, as explained before. Then:
%
\begin{enumerate}

  \item Compute $\mathbf{F}^{n}$ on fluid points, using finite difference approximations.

  \item IBM: On fluid points required for interpolation (i.e. stencil points $e=1,2$ used to define $ex$ interpolation in figure~\ref{Fig:IBCart}) compute the approximation to velocity at time $t_{n+1}$, as defined in the velocity interpolation section. Obtain by interpolation target velocities $\hat{u}_d^{n+1},\hat{v}_d^{n+1}$ on cut-face centroids. Obtain underlying Eulerian face velocities (i.e. $\hat{u}_W^{n+1},\hat{v}_S^{n+1}$ in figure~\ref{Fig:IBCart}) by flux matched averaging, as in equations~\eqref{eq:uE}-\eqref{eq:vN}.

  \item Compute the force $x$ component
  \begin{equation}
    F_x^n = -\frac{\hat{u}_{W}^{n+1}-u_W^n}{\Delta t}  %- \left( \frac{\partial H^{n-1}}{\partial x} \right)_W
  \end{equation}
   at Eulerian locations and assign to vector $\mathbf{F}^{n}$. Same process for $y$ velocity  points.

  \item Solve Poisson equation~\eqref{it:FSPoisson} for $H^{n}$ on Eulerian grid composed of regular and cut-cell underlying Eulerian cells.
  \item Do projection to final velocities~\eqref{it:FSProject}. These velocities respect discretely the divergence constraint on \texttt{GASPHASE} regular and Eulerian cut-cells. Finally, recompute final velocities at cut-face centroids using flux matched averaging.
\end{enumerate}
%





\subsubsection{\&GEOM IB forcing implementation notes}


\subsection{Energy equation, thermodynamic divergence}

\subsubsection{Divergence constraint}

Starting from the sensible enthalpy evolution equation, the divergence of the velocity field can be factored as:
%
\begin{eqnarray}
    ( \nabla \cdot \mathbf{u} )^{th} &=&
    \left[ \frac{1}{\rho c_p T} - \frac{1}{\bar{p}} \right]
    \frac{\partial \bar{p}}{\partial t} + \frac{w \rho_0 g_z}{\rho c_p T} \nonumber \\
    &+& \frac{1}{\rho c_p T} \left[ \dot{q}''' - \nabla \cdot \dot{\mathbf{q}}'' - \mathbf{u} \cdot \nabla (\rho h_s) \right] \nonumber \\
    &+& \frac{1}{\rho} \sum_\alpha \left( \frac{\overline{W}}{W_\alpha} - \frac{h_{s,\alpha}}{c_p T} \right) \left[ \dot{m}_\alpha''' - \nabla \cdot \mathbf{J}_\alpha - \mathbf{u} \cdot \nabla (\rho Y_\alpha) \right] \label{eq:divth}
\end{eqnarray}
%
we call this divergence expression the thermodynamic divergence $( \nabla \cdot \mathbf{u} )^{th}$. The projection scheme for velocities enforces on each cell of the spatial discretization this final divergence on the discrete velocity field.
Next, we look at the finite volume discretization of the non-conservative advection terms involving the velocity $\mathbf{u}$, in the right hand side of equation~\eqref{eq:divth}.

We will use the finite volume version of the previous equation to obtain the volume integrated target thermodynamic divergence on each cut cell volume, and to recompute the thermodynamic divergence on regular gas phase cells adjacent to cut-cells.


\subsubsection{FV discretization of scalar advection terms $\mathbf{u} \cdot \nabla (\phi)$}

Consider the the equality
%
\begin{equation}
    \nabla \cdot \mathbf{\phi u} = \mathbf{u} \cdot \nabla (\phi) + \phi \nabla \cdot \mathbf{u} \label{eq:advforms}
\end{equation}
%
Integrating over the volume of cell $ii$
%
\begin{equation}
    \int_{\Omega_{ii}} {\mathbf{u} \cdot \nabla(\phi)} d\Omega =
    \int_{\Omega_{ii}} \nabla \cdot ({\phi} \mathbf{u}) d\Omega -
    \int_{\Omega_{ii}} \phi \nabla \cdot \mathbf{u} d\Omega \label{eq:flxlim1} \\
\end{equation}
Assuming a flux limited interpolation of the scalar $\phi$ at the cell boundaries, the discrete counterpart is
%
\begin{equation}
    \overline{\mathbf{u} \cdot \nabla(\phi)} \; V_{ii} =
    \sum_{k=1}^{nf_c} (\overline{\phi} \mathbf{u})_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k} \: A_k -
    [\phi]_{ii} \sum_{k=1}^{nf_c} \mathbf{u}_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k} \: A_k \label{eq:flxlim2}
\end{equation}
%
Here the over line states that the cell-centered scalar quantity $\phi \in \{\rho h_s, \rho Y_\alpha\}$ has been interpolated to the cells faces using a limited interpolation scheme. The over line in $\overline{\mathbf{u} \cdot \nabla(\phi)}$, states that this term is computed consistently with the flux limited interpolation adopted. Also, $V_{ii}$ is the volume of cell $ii$, and $A_k$ is the area of a face $k=1,...,nf_c$ of the given cell. This last equation is the finite volume counterpart of equation (B.12) of the FDS technical reference guide.

\subsubsection{Discrete thermodynamic divergence expression}

For cut-cell $ii$ the corresponding discrete volume integrated expression is:
%
\begin{eqnarray}
    ( \nabla \cdot \mathbf{u} )_{ii}^{th} \; V_{ii} &=&
    \left[ \frac{1}{(\rho c_p T)_{ii}} - \frac{1}{\bar{p}_{ii}} \right]
    \frac{\partial \bar{p}_{ii}}{\partial t} \; V_{ii} +
    \frac{w_{ii} \: \rho_0 g_z}{(\rho c_p T)_{ii}} V_{ii} \nonumber \\
    &+& \frac{1}{(\rho c_p T)_{ii}} \left[ \dot{q}''' V_{ii} -
    \sum_{k=1}^{nf_c} \dot{\mathbf{q}}''_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k} \: A_k
    - \overline{\mathbf{u} \cdot \nabla (\rho h_s)} V_{ii} \right] \nonumber \\
    &+& \frac{1}{\rho_{ii}} \sum_\alpha \left( \frac{\overline{W}}{W_\alpha} - \frac{h_{s,\alpha}}{c_p T} \right)_{ii} \left[ \dot{m}_\alpha''' V_{ii} -
    \sum_{k=1}^{nf_c} \mathbf{J}_{\alpha,ii,k} \cdot \hat{\mathbf{n}}_{ii,k} \: A_k
    - \overline{\mathbf{u} \cdot \nabla (\rho Y_\alpha)} V_{ii} \right] \label{eq:divth2}
\end{eqnarray}
%
where the over line terms refer to flux limited interpolation of corresponding scalars, terms defined with subscript $ii$ refer to cell defined quantities, and the scalar diffusive flux $\mathbf{J}_\alpha=- \rho D_\alpha \nabla Y_\alpha$. Also, the vertical velocity $w_{ii}$ \textit{has been interpolated} to the cut-cell centroid.

\subsubsection{Computing diffusive heat flux in FV form:}

Given the species diffusion model adopted the diffusive heat flux vector field is $\dot{\mathbf{q}}_{d,\alpha}''=-h_s \rho D_\alpha \nabla(Y_\alpha)$. The integral of its divergence over a cell $ii$ is:
\begin{equation}
    \int_{\Omega_{ii}} {\nabla \cdot \left(-h_{s,\alpha} \rho D_\alpha \nabla(Y_\alpha) \right)} d\Omega = \sum_{k=1}^{nf_c} \left(-h_{s,\alpha} \rho D_\alpha \nabla(Y_\alpha) \right)_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k}  \: A_k
\end{equation}
where $\left(-h_s \rho D_\alpha \nabla(Y_\alpha) \right)_{ii,k}$ is the mean heat flux defined on boundary surface $k$ of cell $ii$.We employ the same diffusive flux definition as was done on the implicit discretization for scalar transport. Therefore:
\begin{equation}
    \sum_{k=1}^{nf_c} \left(-h_{s,\alpha} \rho D_\alpha \nabla(Y_\alpha) \right)_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k}  \: A_k =
    \sum_{k=1}^{nf_c} -h_{s,\alpha} \left(D_\alpha \nabla(\rho Y_\alpha) -
    \frac{D_\alpha}{\rho} \nabla \rho \; (\rho Y_\alpha) \right)_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k}  \: A_k
\end{equation}


\subsubsection{Computing heat conduction flux in FV form:}


The heat conduction flux vector of the mixture is defined as $\dot{\mathbf{q}}_{kT}''= -k \nabla(T)$, where $k$ is the mixture conductivity. The integral of its divergence over a cell $ii$ is computed similarly as in the previous section:
\begin{equation}
    \int_{\Omega_{ii}} {\nabla \cdot \left(-k \nabla(T)\right)} d\Omega = \sum_{k=1}^{nf_c} \left(-k \nabla(T) \right)_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k}  \: A_k
\end{equation}
where $\left(-k \nabla(T) \right)_{ii,k}$ is the mean heat flux defined on boundary surface $k$ of cell $ii$.





\subsubsection{Cut-cell divergence implementation notes}





\section{Verification Tests}

This section describes some of the test cases used to verify that FDS and Smokeview are computing
and visualizing unstructured geometric objects correctly.  This section will be used as source material for adding
more details to the various guides on how to setup unstructured geometries for use with FDS. (i.e. this section
is not intended to be included verbatim in the FDS User Guide).

Some subsections could become part of the FDS Validation Guide.

\newcommand{\figheightD}{2.75in}

\subsection{Verification of setup checks on broken geometries}

This subsection describes some of the test cases used to verify that FDS correctly enforces on geometries
the quality rules described in Section \ref{triangulated_surfaces_quality}.

The following verification cases contain specific broken geometries in an otherwise correctly
formatted FDS input file:

\begin{itemize}
  \item {\ct geom\_bad\_inconsistent\_normals.fds}: a cube with only one face normal inverted;
  \item {\ct geom\_bad\_inverted\_normals.fds}: a cube with all face normals inverted;
  \item {\ct geom\_bad\_non\_manifold\_edge.fds}: a cube and a box sharing the same non manifold edge;
  \item {\ct geom\_bad\_non\_manifold\_vert.fds}: two cubes sharing the same non manifold vertex; % FIXME currently not enforced
  \item {\ct geom\_bad\_open\_surface.fds}: a cube missing one face, that is an open surface;
\end{itemize}

To allow Firebot to consider the test successful if it catches the error condition,
each of the above input files contain the following parameter in the {\ct MISC}\ line:

{\ct POSITIVE\_ERROR\_TEST = .TRUE.}\

The effect of this parameter is replacing the {\ct ERROR}\ label of the message with a {\ct SUCCESS}\ label.

After Firebot is run, the {\ct geom\_positive\_errors.m}\ Matlab script checks the existence
of the planned error condition in the {\ct CHID.err}\ file, reporting an issue if the condition
is not met.

\subsection{geom\_simple.fds}
Figure \ref{fig:geom_simple} was created using the following GEOM namelist.
The {\ct VERTS}\ keyword is used to specify 3 vertices and the {\ct FACES}\ keyword
is used to specify the indices of the triangular face.

{\small
\begin{verbatim}
&GEOM ID='geom1',
VERTS=0.0,0.0,0.0,1.0,0.5,0.0,1.0,1.0,1.0,
FACES=1,2,3,
SURF_ID='surf1'/
\end{verbatim}
}

\begin{figure}
\begin{center}
\begin{tabular}{c}
 \includegraphics[width=4.0in]{SCRIPT_FIGURES/geom_simple}
  \end{tabular}
\end{center}
 \caption{A triangle created using 3 vertices and 1 face. Case: geom\_simple}
\label{fig:geom_simple}
\end{figure}

%\subsection{geom\_azim.fds}
%Figure \ref{fig:geom_azim} was created using the following GEOM namelists.
%The {\ct AZIM}\ keyword is used to rotate the object about a vertical axis
%centered at $(0,0,0)$.
%
%{\small
%\begin{verbatim}
%&GEOM ID='geom1',
%VERTS=0.0,0.0,0.0, 0.0,1.0,0.0, 0.0,1.0,1.0,
%FACES=1,2,3,AZIM=0.0,SURF_ID='surf1'/
%
%&GEOM ID='geom2',
%VERTS=0.0,0.0,0.0, 0.0,1.0,0.0, 0.0,1.0,1.0,
%FACES=1,2,3,AZIM=45.0,SURF_ID='surf2'/
%
%&GEOM ID='geom3',
%VERTS=0.0,0.0,0.0, 0.0,1.0,0.0, 0.0,1.0,1.0,
%FACES=1,2,3,AZIM=90.0,SURF_ID='surf3'/
%\end{verbatim}
%}
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_azim}
%\end{tabular}
%\end{center}
%\caption{Three triangles generated at different azimuthal angles using the {\ct AZIM}\ keyword. Case: geom\_azim.fds}
%\label{fig:geom_azim}
%\end{figure}
%
%\subsection{geom\_elev.fds}
%Figure \ref{fig:geom_elev} was created using the following GEOM namelists.
%The {\ct ELEV}\ keyword is used to rotate the object with respect to a horizontal plane
%passing through $(0,0,0)$.
%
%{\small
%\begin{verbatim}
%&GEOM ID='geom1',
%VERTS=0.0,0.0,1.0,1.0,0.0,1.0,1.0,1.0,1.0,
%FACES=1,2,3,ELEV=0,XYZ0=0.0,0.0,1.0,SURF_ID='surf1'/
%
%&GEOM ID='geom2',
%VERTS=0.0,0.0,1.0,1.0,0.0,1.0,1.0,1.0,1.0,
%FACES=1,2,3,ELEV=45,XYZ0=0.0,0.0,1.0,SURF_ID='surf2'/
%
%&GEOM ID='geom3',
%VERTS=0.0,0.0,1.0,1.0,0.0,1.0,1.0,1.0,1.0,
%FACES=1,2,3,ELEV=90,XYZ0=0.0,0.0,1.0,SURF_ID='surf3'/
%\end{verbatim}
%}
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_elev}
%\end{tabular}
%\end{center}
%\caption{Three triangles generated at different elevation angles using the {\ct ELEV}\ keyword. Case: geom\_elev.fds}
%\label{fig:geom_elev}
%\end{figure}
%
%\subsection{geom\_scale.fds}
%Figure \ref{fig:geom_scale} was created using the following GEOM namelists.
%The {\ct SCALE}\ keyword is used to change the size of the object.
%
%{\small
%\begin{verbatim}
%&GEOM ID='geom4',
%VERTS=-2.0,0.0,0.0,-1.0,0.0,0.0,-1.0,0.0,1.0,
%FACES=1,2,3,XYZ0=-2.0,0.0,0.0,SCALE=1.0,1.0,1.0,SURF_ID='surf1'/
%
%&GEOM ID='geom4',
%VERTS= 0.0,0.0,0.0,1.0,0.0,0.0,1.0,0.0,1.0,
%FACES=1,2,3,XYZ0= 0.0,0.0,0.0,SCALE=1.0,1.0,2.0,SURF_ID='surf2'/
%
%&GEOM ID='geom4',
%VERTS= 2.0,0.0,0.0,3.0,0.0,0.0,3.0,0.0,1.0,
%FACES=1,2,3,XYZ0= 2.0,0.0,0.0,SCALE=2.0,1.0,1.0,SURF_ID='surf3'/
%\end{verbatim}
%}
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=6.0in]{SCRIPT_FIGURES/geom_scale}
%\end{tabular}
%\end{center}
%\caption{Three triangles generated at different scale sizes using the {\ct SCALE}\ keyword. Case: geom\_scale.fds}
%\label{fig:geom_scale}
%\end{figure}

\subsection{geom\_obst.fds}
Figure \ref{fig:geom_obst} was created using the following GEOM namelist.
The {\ct XB}\ keyword is used the same way to specify a block as on
a {\ct OBST}\ or {\ct VENT}\ line.

{\scriptsize
\begin{verbatim}
&GEOM ID='obst',SURF_ID='surf1', XB=-0.6,0.6,-0.6,0.6,-0.2,0.2,AZIM=90.0,ELEV=30.0 /
\end{verbatim}
}

\begin{figure}
\begin{center}
\begin{tabular}{c}
 \includegraphics[height=\figheightD]{SCRIPT_FIGURES/geom_obst}
  \end{tabular}
\end{center}
 \caption{A block generated using the {\ct XB}\ keyword.  The block is refined automatically to be consistent with the underlying grid resolution. Case: geom\_obst.fds}
\label{fig:geom_obst}
\end{figure}

\subsection{geom\_sphere1a.fds,...,geom\_sphere3f.fds}
Figure \ref{fig:geom_sphere} was created using {\ct LEVEL=n}\
in the following GEOM namelist
where {\ct n}\ ranges from 0 to 5.

{\scriptsize
\begin{verbatim}
&GEOM ID='sphere',SURF_ID='surf1',SPHERE_RADIUS=0.5,N_LEVELS=n,SPHERE_ORIGIN=0.0,0.0,0.0 /
\end{verbatim}
}

The {\ct N\_LEVELS}\ keyword is used
to specify the resolution of the sphere, larger {\ct N\_LEVELS}\ values result
in a more highly resolved sphere.{\ct LEVEL=0} produces a 20 sided sphere approximation, an icosahedron .
{\ct LEVEL=n} produces a sphere approximation with four times as many triangles as the
sphere produced with {\ct LEVEL=n-1}.
The {\ct SPHERE\_RADIUS}\ and {\ct SPHERE\_ORIGIN}\ keywords are used to specify
the size and location of the sphere.  This discretization technique results in equilateral triangles at each recursion level {(\em check to make sure this statement is true - when verified remove this comment)}.

\begin{figure}
\begin{center}
\begin{tabular}{cc}
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere1a}&
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere1b}\\
 LEVEL=0&LEVEL=1\\
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere1c}&
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere1d}\\
 LEVEL=2&LEVEL=3\\
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere1e}&
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere1f}\\
 LEVEL=4&LEVEL=5\\
  \end{tabular}
\end{center}
 \caption{Recursive sphere discretizations.  A sphere at a given level is
 obtained by splitting each triangle from the previous level into four parts and renormalizing added vertices. Cases: geom\_sphere1a.fds,...,geom\_sphere1f.fds}
\label{fig:geom_sphere}
\end{figure}

Figure \ref{fig:geom_sphere2a} was created by using {\ct N\_LAT}\ and {\ct N\_LONG}\ keywords in a {\ct GEOM}\ namelist
to split a sphere in latitudinal (north/south) and longitudinal (east/west) directions respectively. The minimum values of {\ct N\_LAT}\ and {\ct N\_LONG} permitted are 3 and 6.  This discretization technique results in triangles with high aspect ratios near the poles at higher discretization levels.

\begin{figure}
\begin{center}
\begin{tabular}{cc}
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere3a}&
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere3b}\\
 N\_LAT=3,N\_LONG=6&N\_LAT=6,N\_LONG=12\\
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere3c}&
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere3d}\\
 N\_LAT=12,N\_LONG=24&N\_LAT=24,N\_LONG=48\\
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere3e}&
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_sphere3f}\\
 N\_LAT=48,N\_LONG=96&N\_LAT=96,N\_LONG=192\\
  \end{tabular}
\end{center}
 \caption{Latitude/Longitude sphere discretizations.  Spheres are
 split in longitudinal (east/west) and latitudinal (north/south) directions using the {\ct N\_LAT}\ and {\ct N\_LONG} keywords. Cases: geom\_sphere3a.fds,...,geom\_sphere3f.fds}
\label{fig:geom_sphere2a}
\end{figure}

%\geominput{geom_terrain.fds}
\subsection{geom\_terrain.fds}
Figure \ref{fig:geom_terrain} was created using the {\ct ZVALS}\ keyword.

\begin{figure}
\begin{center}
\begin{tabular}{c}
 \includegraphics[height=\figheightD]{SCRIPT_FIGURES/geom_terrain}
  \end{tabular}
\end{center}
 \caption{Elevations are defined on a rectangular array of grid points using the {\ct ZVALUES}\ keyword.  Case: gridgeom\_terrain.fds}
\label{fig:geom_terrain}
\end{figure}

\subsection{geom\_texture.fds}
Figure \ref{fig:geom_texture} was created using the following SURF and GEOM namelists.
The {\ct TEXTURE\_MAP}\ keyword on the {\ct SURF}\ is used to specify the name of the image
file used to texture map the geometric object.

{\small
\begin{verbatim}
&SURF ID='surf1',TEXTURE_MAP='nistleft.jpg',TEXTURE_WIDTH=0.6,TEXTURE_HEIGHT=0.2,COLOR='BLUE' /

&GEOM ID='texture',
VERTS=0.0,0.0,0.0, 1.0,0.0,0.0, 1.0,1.0,0.0,
FACES=1,2,3,SURF_ID='surf1'/
\end{verbatim}
}

\begin{figure}
\begin{center}
\begin{tabular}{c}
 \includegraphics[width=4.0in]{SCRIPT_FIGURES/geom_texture}
  \end{tabular}
\end{center}
 \caption{A texture is applied to a triangle. Case: geom\_texture.fds}
\label{fig:geom_texture}
\end{figure}

\subsection{geom\_texture2.fds}
Figure \ref{fig:geom_texture2} was created using the following SURF and GEOM namelists.
The {\ct TEXTURE\_MAP}\ keyword on the {\ct SURF}\ is used to specify the name of the image
file used to texture map the geometric object. This example has two textures.

{\small
\begin{verbatim}
&SURF ID='surf1' TEXTURE_MAP='nistleft.jpg',TEXTURE_WIDTH=0.6,TEXTURE_HEIGHT=0.2,COLOR='BLUE' /
&SURF ID='surf2',TEXTURE_MAP='grass.jpg',TEXTURE_WIDTH=0.6,TEXTURE_HEIGHT=0.2,COLOR='GREEN' /

&GEOM ID='texture',
VERTS=0.0,0.0,0.0, 1.0,0.0,0.0, 1.0,1.0,0.0,
FACES=1,2,3,SURF_ID='surf1'/

&GEOM ID='texture2',
VERTS=0.0,0.0,0.0, 1.0,1.0,0.0, 0.0,1.0,0.0,
FACES=1,2,3,SURF_ID='surf2'/
\end{verbatim}
}

\begin{figure}
\begin{center}
\begin{tabular}{c}
 \includegraphics[width=4.0in]{SCRIPT_FIGURES/geom_texture2}
  \end{tabular}
\end{center}
 \caption{Two texture maps are applied to two separate triangles.  Case: geom\_texture2.fds}
\label{fig:geom_texture2}
\end{figure}

\subsection{geom\_texture3a.fds, geom\_texture3b.fds, geom\_texture4a.fds, geom\_texture4b.fds}
The spheres in Figure \ref{fig:geom_texture3} were created using SURF and GEOM namelists similar to the following.

{\small
\begin{verbatim}
&SURF ID='surf1',TEXTURE_MAP='sphere_cover_03.png'COLOR='BLUE',
      TEXTURE_WIDTH=1.0,TEXTURE_HEIGHT=1.0/

&GEOM ID='sphere1',SURF_ID='surf1',
      N_LAT=50,N_LONG=50,SPHERE_RADIUS=0.25,SPHERE_ORIGIN=0.25,0.5,0.5,
      TEXTURE_ORIGIN=0.25,0.5,0.5,TEXTURE_MAPPING='SPHERICAL' /

&GEOM ID='sphere2',SURF_ID='surf3',
      N_LEVELS=3,SPHERE_RADIUS=0.25,SPHERE_ORIGIN=0.75,0.5,0.5,
      TEXTURE_ORIGIN=0.75,0.5,0.5,TEXTURE_MAPPING='SPHERICAL' /
\end{verbatim}
}

The {\ct TEXTURE\_MAP}\ keyword is used to specify the name of the image
file applied to the geometric object. The {\ct SPHERICAL}\ setting indicates that that the texture map image
is applied to the object using spherical coordinates.



\begin{figure}
\begin{center}
\begin{tabular}{cc}
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_texture3a}&
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_texture3b}\\
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_texture4a}&
 \includegraphics[width=3.5in]{SCRIPT_FIGURES/geom_texture4b}\\
 lat/long discretization&recursive discretization
  \end{tabular}
\end{center}
 \caption{Texture maps applied to spheres discretized using two different methods.
 The two sphere on the left are discretized by splitting the sphere along longitudinal (east/west) and latitudinal (north/south) directions.
 The two spheres on the right are discretized recursively starting with an icosahedron ( 20 sided polyhedron).  Texture maps applied to the top two sphere have 8 segments. Texture map applied to the bottom two spheres have 5 segments.  Cases: geom\_texture3a.fds, geom\_texture3b.fds, geom\_texture4a.fds, geom\_texture4b.fds}
\label{fig:geom_texture3}
\end{figure}

\subsection{geom\_arch.fds}
Figure \ref{fig:geom_arch} was created using VERTS and FACES keywords.

\begin{figure}
\begin{center}
\begin{tabular}{c}
 \includegraphics[width=4.0in]{SCRIPT_FIGURES/geom_arch}
  \end{tabular}
\end{center}
 \caption{An arch structure created with FACES and VERTS keywords.
 Case: geom\_arch.fds . }
\label{fig:geom_arch}
\end{figure}

%\subsection{geom\_time.fds}
%Figure \ref{fig:geom_time} was created using the following GEOM namelist.
%The {\ct AZIM}\ keyword is used to specify the initial rotational angle
%and the {\ct AXIM\_DOT}\ is used to specify that rate at which
%the azimuthal angle changes, in this  case 0.36 deg/s.
%
%{\small
%\begin{verbatim}
%&GEOM ID='geom4',
%      VERTS=0.0,0.0,0.0,0.0,1.0,0.0,0.0,1.0,1.0,
%      FACES=1,2,3,AZIM=0.0,AZIM_DOT=0.36,SURF_ID='surf1'/
%\end{verbatim}
%}
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{ccc}
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time_050}&
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time_100}&
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time_150}\\
% \SI{50.0}{s}&\SI{100.0}{s}&\SI{150.0}{s}
%  \end{tabular}
%\end{center}
% \caption{A triangle rotating about the scene center. Case: geom\_time.fds}
%\label{fig:geom_time}
%\end{figure}
%
%\subsection{geom\_time2.fds, geom\_time3.fds, geom\_time4.fds}
%Figure \ref{fig:geom_time2} was created using the following GEOM namelists.
%The {\ct AZIM}\ keyword is used to specify the initial rotational angle
%and the {\ct AXIM\_DOT}\ is used to specify that rate at which
%the azimuthal angle changes, in this  case 0.72 deg/s.
%The {\ct blade}\ object is used as components to create four blade objects
%{\ct blade1}, {\ct blade2}, {\ct blade3} and {\ct blade4} of the
%{\ct prop}\  object. The blades are positioned within the {\ct prop} using the {\ct DAZIM}\ keyword. Figures \ref{fig:geom_time3} and \ref{fig:geom_time4} were created similarly using rotation keywords to incline the geometric objects.
%
%{\small
%\begin{verbatim}
%&GEOM ID='blade',AZIM=0.0,AZIM_DOT=0.72
%VERTS=0.0,0.0,0.0,1.0,0.0,0.0,1.0,0.0,1.0,
%FACES=1,2,3,COMPONENT_ONLY=.TRUE.,SURF_ID='surf4'/
%
%&GEOM ID='blade1',GEOM_IDS(1)='blade',SURF_ID='surf1',COMPONENT_ONLY=.TRUE. /
%&GEOM ID='blade2',GEOM_IDS(1)='blade',SURF_ID='surf2',COMPONENT_ONLY=.TRUE. /
%&GEOM ID='blade3',GEOM_IDS(1)='blade',SURF_ID='surf3',COMPONENT_ONLY=.TRUE. /
%&GEOM ID='prop' GEOM_IDS(1)='blade1',DAZIM(1)=0.0,
%GEOM_IDS(2)='blade2',DAZIM(2)=90.0,
%GEOM_IDS(3)='blade1',DAZIM(3)=180.0,
%GEOM_IDS(4)='blade3',DAZIM(4)=270.0,
%\end{verbatim}
%}
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{ccc}
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time2_050}&
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time2_100}&
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time2_150}\\
% \SI{50.0}{s}&\SI{100.0}{s}&\SI{150.0}{s}
%  \end{tabular}
%\end{center}
% \caption{Four triangles defined as a group rotating about the scene center. A Case: geom\_time2.fds}
%\label{fig:geom_time2}
%\end{figure}
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{ccc}
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time3_050}&
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time3_100}&
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time3_150}\\
% \SI{50.0}{s}&\SI{100.0}{s}&\SI{150.0}{s}
%  \end{tabular}
%\end{center}
% \caption{Four inclined triangles defined as a group rotating and about the scene center. A Case: geom\_time3.fds}
%\label{fig:geom_time3}
%\end{figure}
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{ccc}
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time4_050}&
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time4_100}&
% \includegraphics[width=2.0in]{SCRIPT_FIGURES/geom_time4_150}\\
% \SI{50.0}{s}&\SI{100.0}{s}&\SI{150.0}{s}
%  \end{tabular}
%\end{center}
% \caption{Two sets of four triangles rotating about the scene center. Each set of 4 triangles are defined as a group. Case: geom\_time4.fds}
%\label{fig:geom_time4}
%\end{figure}

%\subsection{geom\_slice.fds}
%The slice file in Figure \ref{fig:geom_slice} was generated using the {\tt SLICETYPE}\ keyword on a {\tt SLCF}\ namelist statement as in
%{\footnotesize
%\begin{verbatim}
%&SLCF PBY=0.75,QUANTITY='TEMPERATURE',CELL_CENTERED=.TRUE., SLICETYPE='EXCLUDEGEOM' /.
%\end{verbatim}
%}
%
%\noindent The {\tt SLICETYPE}\ keyword is used to specify the type of slice generated.  For example, a slice consisting of only cutcells, a slice  including the effects of immersed geometry, a slice ignoring immersed geometry and a slice output using the original slice file format.  These choices are documented in Table \ref{table:slicetype}.
%
%\begin{table}
%\caption{Description of {\tt SLICETYPE}\ values used on a {\tt SLCF}\ namelist statement }
%\begin{tabular}{|l|l|}
%  \hline
%Keyword& Description\\
%\hline
%STRUCTURED&generate slice using original  file format\\
%CUTCELL&output only triangle faces associated with  cutcells\\
%INCLUDEGEOM&include immersed geometry when generating a slice file\\
%EXCLUDEGEOM&ignore immersed geometry when generating a slice file\\
%\hline
%\end{tabular}
%\label{table:slicetype}
%\end{table}
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{cc}
% \includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_slice_slice_000}&
% \includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_slice_slice_020}\\
% \SI{0.0}{s}&\SI{20.0}{s}
%  \end{tabular}
%\end{center}
% \caption{Slice file generated using the geometry file format Case: geom\_slice.fds.  The two triangles are an immersed object with boundary file temperature data displayed on its surface.}
%\label{fig:geom_slice}
%\end{figure}

%\subsection{geom\_sphere2.fds}
%Figure \ref{fig:geom_sphere2b} tests motion of a scaled,  moving, rotating sphere. The sphere was scaled using {\ct SCALE=0.6,0.4,0.2}\, placed initially at the right side of the domain using {\ct XYZ=10.0,0.5,1.0}, translated to the left and down using {\ct XYZ\_DOT=-0.01,0.0,-0.001} and rotated using {\ct GROTATE\_DOT=0.36}.  The axis of rotation was specified using {\ct GAXIS=0.0,1.0,0.0}.
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=6.0in]{SCRIPT_FIGURES/geom_sphere2_000}\\
%\includegraphics[width=6.0in]{SCRIPT_FIGURES/geom_sphere2_200}\\
%\includegraphics[width=6.0in]{SCRIPT_FIGURES/geom_sphere2_400}\\
%\includegraphics[width=6.0in]{SCRIPT_FIGURES/geom_sphere2_600}\\
%\includegraphics[width=6.0in]{SCRIPT_FIGURES/geom_sphere2_800}\\
%\includegraphics[width=6.0in]{SCRIPT_FIGURES/geom_sphere2_1000}
%  \end{tabular}
%\end{center}
% \caption{Motion of a scaled, moving, rotating sphere.  Case: geom\_sphere2.fds}
%\label{fig:geom_sphere2b}
%\end{figure}

%\subsection{geom\_sphere\_fire.fds}
%Figure \ref{fig:geom_sphere_fire} tests visualization of a boundary file data on a sphere.

%\begin{figure}
%\begin{center}
%\begin{tabular}{cc}
%\includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere_fire_000}&
%\includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere_fire_020}\\
%\SI{0.0}{s}&\SI{2.0}{s}\\
%\includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere_fire_040}&
%\includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere_fire_060}\\
%\SI{4.0}{s}&\SI{6.0}{s}\\
%\includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere_fire_080}&
%\includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere_fire_100}\\
%\SI{8.0}{s}&\SI{10.0}{s}\\
% \end{tabular}
%\end{center}
%\caption{Boundary file data visualized on a sphere.  Case: geom\_sphere\_fire.fds}
%\label{fig:geom_sphere_fire}
%\end{figure}

%Figure \ref{fig:geom_sphere2_fire} tests visualization of a boundary file data on a moving sphere.

%\begin{figure}
%\begin{center}
%\begin{tabular}{cc}
% \includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere2_fire_000}&
% \includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere2_fire_020}\\
% \SI{0.0}{s}&\SI{2.0}{s}\\
% \includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere2_fire_040}&
% \includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere2_fire_060}\\
% \SI{4.0}{s}&\SI{6.0}{s}\\
% \includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere2_fire_080}&
% \includegraphics[width=3.0in]{SCRIPT_FIGURES/geom_sphere2_fire_100}\\
% \SI{8.0}{s}&\SI{10.0}{s}\\
%  \end{tabular}
%\end{center}
% \caption{Boundary file data visualized on a moving sphere.  Case: geom\_sphere2\_fire.fds}
%\label{fig:geom_sphere2_fire}
%\end{figure}

\clearpage

\subsection{Cut-cell definition around immersed cone: split cells, piercing and alignment}

\label{cone_test}
We test the computational geometry algorithm on its treatment of degenerate cases~\cite{Berger:2017}, namely, split cells within a cartesian cell, cartesian cell piercing, and body surface-grid alignment. These cases are expected to be found in general. A cone of base diameter $D_B$ and height $H=4D_B$ is placed centered in a Cartesian domain of size $x=[-1.5D_B:1.5D_B]$, $y=[-1.5D_B:1.5D_B]$, $z=[-2D_B:6D_B]$, as shown in figure~\ref{Fig:cone_cut_cells}a. The cone is composed of 316 vertices and 628 triangles. A single Cartesian mesh is used with $n_x=5$, $n_y=5$ and $n_z=12$ cells in each corresponding direction. 

A center column of cartesian cells and corresponding cut-cells is shown in figure~\ref{Fig:cone_cut_cells}b. The proposed method correctly defines cut-faces and polyhedrons in three situations, namely: Solid piercing a cartesian cell, split cut-cells, and cells defined from the cone base aligned with grid plane $z=0$.


 %Figure:
\begin{figure}[h]
      \centering
      \includegraphics[trim = 80mm 10mm 80mm 10mm, clip,
      width=0.325\linewidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/cone_1mesh_a.png} 
      \includegraphics[trim = 80mm 10mm 80mm 10mm, clip,
      width=0.325\linewidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/cone_1mesh_b.png} 
      \put(-310,240){(a)}    
      \put(-120,240){(b)}
      \put(-30,185){$\leftarrow$  pierced cell}    
      \put(-30,90){$\leftarrow$  split cells}
      \put(-30,40){$\leftarrow$ aligned face}
      \caption{Polyhedra computation around immersed cone: (a) problem setup and grid on a $y$ plane. (b) Cut-cells on mid cell column along $z$ direction: \textit{(pierced cell)} cut-cell formed by cartesian cell pierced by object; \textit{(split cells)} cartesian cell complement of body defines 4 split cells; and \textit{(aligned face)} cut-cell of same size of regular cell contains boundary cut-faces aligned with $z=0$ grid plane.}
	\label{Fig:cone_cut_cells}
\end{figure}


Three integral conservation tests are defined for a completely immersed geometry problem. These are namely, relative difference of areas computed from the surface triangulation or from the computed boundary cut-faces $|A_{B,\tau}-A_{B,CF}|/A_{B,\tau}$; relative difference of $G$ region volumes computed from the body triangulation and from the cut-cells and $G$ region regular cells $|V_{GB}-V_{GRC}|/V_{GB}$; and finally, difference on the $G$ region centroid computed from body centroid or cut-cells and $G$ cells centroids $|\mathbf{x}_{GB}-\mathbf{x}_{GRC}|_\infty$. See Appendix~\ref{info:constests} for a formal definition of these tests.

The results of the three tests on this cone problem can be seen on table~\ref{tab:cone_res}. We note that cone area and $G$ region volume computed either through the solid triangulation or through the polyhedra produced by the method match up to relative values of $10^{-15}$, commensurate with the floating point tolerance $\epsilon=10^{-12}L_z$ used. In practice, this level of correspondence guarantees conservation as far as the volume mesh is concerned. The centroid of the $G$ region computed through the two different avenues matches up to $10^{-16}D_B$, and is an excellent integral test of the proposed method to compute both cut-face and cut-cell centroids. Note that centroids along $x$ and $y$ direction are not zero, due to the asymmetry of the geometry triangulation employed.   

%
\begin{table}[h]
  \begin{tabular}{c|c|c}
    \hline \hline
      $A_{B,\tau}$ & $A_{B,CF}$ & $\frac{|A_{B,\tau}-A_{B,CF} |}{A_{B,\tau}}$ \\ \hline
       \texttt{0.7104E+01} & \texttt{0.7104E+01} & \texttt{0.6252E-15} \\
    \hline \hline
      $V_{GB}$ & $V_{GRC}$ & $\frac{|V_{GB}-V_{GRC}|}{V_{GB}}$ \\  \hline
      \texttt{0.7096E+02} & \texttt{0.7096E+02} & \texttt{0.2212E-14} \\
    \hline \hline
    $\mathbf{x}_{GB}$ & $\mathbf{x}_{GRC}$ & $|\mathbf{x}_{GB}-\mathbf{x}_{GRC}|_\infty$ \\
    \hline \hline
    \texttt{(-0.306E-6,0.586E-6,0.202E+1)} &
    \texttt{(-0.306E-6,0.586E-6,0.202E+1)} &
    \texttt{0.1856E-15} \\
    \hline
  \end{tabular}
  \caption{Computation results for conservation tests on polyhedrons defined on immersed cone problem.} 
  \label{tab:cone_res}
\end{table}
%


%\subsection{Cut-cell definition around sphere: serial performance and parallel scaling}


\subsection{Explicit-implicit time integration for scalar transport}


\subsubsection{Temporal Error analysis for variable density projection}

\label{sec:saad_cc_temporal_error}

This case is a variation of problem \texttt{sec:saad\_temporal\_error} for a domain that includes an imposed cut-cell region, where the cut-cells have same geometry as regular cells. The problem allows us to verify the temporal accuracy of our scalar transport implementation for cut-cells. The problem is defined by the advection of a non-reacting mixture of two fluids with different densities (10:1 ratio) in one cartesian direction. Zero diffusivity is prescribed. For the problem definition and temporal error computation details, we refer the reader to problem \texttt{sec:saad\_temporal\_error}.
The domain is defined by a mesh containing $N_x=512$ cells in the $x$ direction, where the domain $x=[-1,1]$. Also, four cells are defined in the $y$ direction as required by the FFT Poisson solver. Note that the problem defined is one-dimensional in $x$.
On this domain, cut-cells with same shape and size of regular cells are defined in the space $x=[-0.5,0.5]$. The test verifies that the implementation details of this cut-cell region does not affect the temporal order of accuracy seen on the default verification case of section \texttt{sec:saad\_temporal\_error}.

\title{Temporal order of accuracy}

Six cases are defined for this test:
%
\begin{enumerate}
 \item[a]  Three correspond to default FDS DNS flux limiter and default explicit Strong Stability Preserving Runge Kutta (SSPRK2) time integration in cells not covered by cut-cells, plus Godunov flux limiter and explicit SSPRK2 time integration in cut-cell region, for CFL numbers $0.25$, $0.125$, $0.0625$. These cases are defined in {\ct /Complex\_Geometry/saad\_CC\_explicit\_*.fds}.

  \item[b] The other three cases correspond to default flux limiter and time integration in cells not covered by cut-cells, plus Godunov flux limiter and implicit time integration in cut-cell region for the same CFL numbers. The implicit integration on the cut-cell region is done applying an Backward Euler step in the predictor, and Trapezoidal Rule on the corrector step. The imposition of boundary conditions for the implicit region in domain and explicit-implicit (EXIM) boundaries follows the technical notes in section \texttt{sec:exim\_scalar\_transport} of the Technical Reference Guide. These cases are defined in {\ct /Complex\_Geometry/saad\_CC\_implicit\_*.fds}.

\end{enumerate}
%

Same procedure as explained in section  \texttt{sec:saad\_temporal\_error} is used to define the temporal order of accuracy $p$. As noted there, this procedure filters out the spatial error from the computed temporal convergence. Initial and final fields for temperature and mixture fraction for the smallest integration time step are shown in the top figures of Fig.~\ref{fig:saad_cc_temporal_order}. Fields advect to the right with unitary velocity. The two bottom figures show the order $p$ computed point wise in the $x$ direction for explicit and implicit time integration on the cut-cell region. Deviations from $p \simeq 2$ are dependent on the formula used to compute $p$, as explained in section~\texttt{sec:saad\_temporal\_error}.
The $l_2$ norm of $p$ computed pointwise for explicit integration in cut-cells is \input{../FDS_Verification_Guide/SCRIPT_FIGURES/saad_CC_explicit_l2_norm}\!, and for implicit time integration of scalars on these is  \input{../FDS_Verification_Guide/SCRIPT_FIGURES/saad_CC_implicit_l2_norm}\!.

% Figure Saad_CC:
\begin{figure}[ht]
\centering
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/saad_CC_explicit_Z.pdf}
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/saad_CC_explicit_rho.pdf}
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/saad_CC_explicit_temporal_order_rho.pdf}
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/saad_CC_implicit_temporal_order_rho.pdf}
\caption[The {\ct saad CC} temporal order test case]{Temporal order for a variable-density projection on domain including cut-cell region.  (Upper-left) Initial and final field for mixture fraction, explicit time integration on cut-cells.  (Upper-right) Initial and final field for density, explicit time integration on cut-cells.  (Lower-left) $p$ computed pointwise for density at final time, explicit integration on cut-cells. (Lower-right) $p$ computed pointwise for density at final time, implicit time integration on cut-cells.  Fluctuations are due to degenerate points in the formula for $p$.}\label{fig:saad_cc_temporal_order}
\end{figure}


\subsubsection{Variable-Density manufactured solution}

\label{sec:shunn_cc_mms}


We set to test the manufactured solution test case~\texttt{sec:shunn\_mms} for scalar transport, including an internal region of the two dimensional domain where cut-cells with same shape and size as regular Cartesian cells are artificially defined. This problem allows us to check the overall accuracy (spatial+temporal) of the cut-cell implementation for scalar transport for both explicit and implicit integration on the cut-cell region.
The parameters of the test follow the ones defined by Shunn et al.~\cite{Shunn:2012} and described in section~\texttt{sec:shunn\_mms}.

The treatment of sources for scalar transport and momentum equations due to the manufactured solution is essentially the same as defined in section~\texttt{sec:shunn\_mms}. Therefore, source terms from manufactured solution are treated following the base SSPRK2 integration method defined in FDS.
The domain is defined on a square of length ($L=2$ m), with periodicity prescribed in its boundaries. A smaller square of length ($L=1$m) centered on the previous is filled with cut-cells of regular cell geometry. See figure~\texttt{fig:smokeview\_fig}.


% Plot of a slice of temperature on the whole domain with the cut-cell edges showing the cut-cell region:
% Figure:


Simulations were performed for $N = \{32, 64, 128, 256, 384\}$, where $N$ is the number of cells in each direction, using an adaptive time step satisfying a CFL number of 0.25. The scalar transport equations for both species defined on section~\texttt{sec:shunn\_mms}, are solved using the GODUNOV flux limiter and centered interpolation. Figure~\ref{fig:shunn_cc_accuracy_order}, presents the $L_2$ error norm defined on the whole domain (regular+cut-cells) at time $t = 0.9$ as a function of grid spacing, $\Delta x = L/N$. Both explicit and implicit time integration in the cut-cell region produce second-order accurate solutions for density, mixture fraction, and velocity, when centered interpolation is used to define the advective mass flux. As expected, these variables exhibit asymptotic first order accuracy when the GODUNOV flux limiter is used. Note that central differences are used on the diffusive terms of these equations, rendering a second order accurate approximation of diffusion on these regular shaped cells. The pressure $H$ is, in same fashion as seen on section~\texttt{sec:shunn\_mms}, first order accurate. This is consistent with the projection scheme used for time integration of velocities used by FDS.


% Figure Shunn_CC:
\begin{figure}[ht]
\centering
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/shunn_cc_exp_cen_mms_convergence.pdf}
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/shunn_cc_exp_gdv_mms_convergence.pdf}
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/shunn_cc_imp_cen_mms_convergence.pdf}
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/shunn_cc_imp_gdv_mms_convergence.pdf}
\caption[The {\ct Shunn CC} accuracy order test case]{Overall accuracy for manufactured solution variable-density projection on domain including cut-cell region. The $L_2$ error defined on the whole domain at time $t = 0.9$ s is plotted as a function of grid spacing for $N=\{32,64,128,256,384\}$ points in each direction. Two flux limiting options for scalar advection are used: centered interpolation and Godunov.  (Upper-left) Convergence for explicit scalar time integration and centered interpolation for advection.  (Upper-right) Convergence for explicit scalar time integration and Godunov flux limiter.  (Lower-left) Convergence for implicit scalar time integration (Backward Euler step in predictor and Trapezoidal rule in corrector step), and  centered interpolation for advection. (Lower-right) Convergence for implicit scalar time integration and Godunov flux limiter.  }\label{fig:shunn_cc_accuracy_order}
\end{figure}





\subsection{Immersed boundary method on manufactured solution around rotated cube}

Consider the solution of the following constant density, incompressible flow, momentum and scalar transport equations:
%
\begin{eqnarray}
  \frac{\partial \rho \mathbf{u}(\mathbf{x},t)}{\partial t} &=& -\left[ \mathbf{F_v}(\mathbf{u},\mathbf{x},t) + \boldsymbol{\nabla} p(\mathbf{x},t) + \mathbf{fe_v}(\mathbf{u},p) \right] \label{eq:mmsMom} \\
  \frac{\partial \rho Z}{\partial t} &=& -\left[   \mathbf{F_z}(\mathbf{u},\mathbf{x},t) +  \mathbf{fe_Z}(\mathbf{u},Z)   \right]
  \label{eq:mmsZ}
\end{eqnarray}
%
which are extended by the manufactured solution (MMS) sources $\mathbf{fe_v}(\mathbf{u},p)$ and $\mathbf{fe_Z}(\mathbf{u},Z)$ respectively. The vectors $\mathbf{F_v}=-\mu \nabla^2(\mathbf{u}) + \nabla \cdot \left( \rho \mathbf{u} \otimes \mathbf{u} \right)$ and $\mathbf{F_z}= \nabla \cdot \left( \rho Z \mathbf{u} - \rho D_z \nabla Z \right) $ contain the corresponding advection and diffusion operators.

Now, consider in figure~\ref{Fig:ManuSoln}, a local set of fixed axes $x',y'$ and a manufactured solution defined on the fluid region $\Omega - \Omega_1$ of the form:
%
\begin{eqnarray}
   u'(\mathbf{x}_{Bp}',t) &=& -\sin(t) \; \sin(kx')^2 \; \sin(2ky') \label{eq:upmf} \\
   v'(\mathbf{x}_{Bp}',t) &=&   \sin(t) \; \sin(2kx') \; \sin(ky')^2 \label{eq:vpmf} \\
   p(\mathbf{x}_{Bp}',t) &=& \frac{\sin(t)}{4} \left( 2+\cos(2kx') \right) \left( 2+\cos(2ky') \right) - \sin(t) \label{eq:pmf} \\
   Z(\mathbf{x}_{Bp}',t) &=& \frac{A_Z \sin(t)}{3} \left[ \left( 1-\cos(2k(x'-\frac{\pi}{2})) \right)
                                                                                   \left( 1-\cos(2k(y'-\frac{\pi}{2})) \right) - 1 \right] + Z_{mean} \label{eq:qmf}
\end{eqnarray}
%
which is similar on velocities and pressure to the functions given in reference~\cite{Weinan:1}.  These functions imply zero initial conditions, a divergence free velocity field and the following
%
\begin{equation}
   u'(\mathbf{x}',t)=v'(\mathbf{x}',t)=0, \; \; \frac{\partial p}{\partial n} = \frac{\partial Z}{\partial n}=0, \; \;
   for \; \; kx'=ky'= m \pi \; ; \; m \in \mathbb{N}
\end{equation}
%
which can be used as \textit{no-slip wall} conditions in square boundaries defined by the integer $m$. Indeed, one could be interested in the flow field inside the square boundaries, or outside.

 The MMS sources $\mathbf{fe_v}$ and $\mathbf{fe_Z}$ are computed analytically as functions of $(\mathbf{x},t)$ replacing definitions~\eqref{eq:upmf}-\eqref{eq:qmf} in the governing equations~\eqref{eq:mmsMom}-\eqref{eq:mmsZ} and computing symbolically the corresponding terms. The Matlab script that performs said derivations can be found in the fds repository in \textbf{DEFINE LOCATION of matlab script}.



We take values $m=0,1$ and wavenumber $k=1$, and consider the external flow to the resulting square immersed body, as seen in figure~\ref{Fig:ManuSoln}. The velocity field function inside the immersed object $\Omega_1$ (i.e. solid region) is assumed zero at all times, consistent with the fixed rigid body hypothesis.
%
% figure
\begin{figure}[h]
      \centering
      \includegraphics[trim = 40mm 20mm 40mm 20mm, clip,
       width=0.6\linewidth]{../../../fig/fds/Geometry_Figures/RotCube_ManuSolnSketch.png}
      \caption{Reference frames for Eulerian ($E$) and immersed body local system ($B$), and vectors relating a generic point $p$ on the fluid region $\Omega-\Omega_1$.}
	\label{Fig:ManuSoln}
\end{figure}
%
An angular parameter $\alpha$ defines the orientation of the local axes $x',y'$ respect to the Eulerian frame coordinates $x,y$. The center of $\Omega_1$ in the Eulerian frame is situated in position $\mathbf{x}_{Ec}=(\pi,\pi)$. The rotation matrix from Eulerian $E$ to local frame $B$ coordinates is
%
\begin{equation}
  [\mathbf{T}_{BE}] = \left[ \begin{array}{c c}
                                          \cos(\alpha) & \sin(\alpha) \\
                                         -\sin(\alpha)  & \cos(\alpha) \end{array} \right]
\end{equation}
%
In the local system, the center of the square in local coordinates is given by $\mathbf{x}_{Bc}'=(\pi/2,\pi/2)$. Then, the location of a general point $p$ in the local system is:
%
\begin{equation}
    \mathbf{x}_{Bp}'= \mathbf{x}_{Bc}' +\mathbf{x}_{cp}'
\end{equation}
%
The same point $p$ location can be defined in terms of the $E$ frame as $\mathbf{x}_{Ep}= \mathbf{x}_{Ec} +\mathbf{x}_{cp}$ (Eulerian coordinates). Then
%
\begin{equation}
    \mathbf{x}_{cp}'=[\mathbf{T}_{BE}] \mathbf{x}_{cp}=[\mathbf{T}_{BE}] \left( \mathbf{x}_{Ep} - \mathbf{x}_{Ec} \right)
\end{equation}
%
and finally
%
\begin{equation}
    \mathbf{x}_{Bp}'= \mathbf{x}_{Bc}' +[\mathbf{T}_{BE}] \left( \mathbf{x}_{Ep} - \mathbf{x}_{Ec} \right)
\end{equation}
%
This last equation is employed to translate a given point $\mathbf{x}_{Ep}$ in the Eulerian frame (with Eulerian frame coordinates $x,y$) to the local frame coordinates where equations~\eqref{eq:upmf}-\eqref{eq:qmf} are defined. This transformation is used to define the initial conditions in the local frame, for points collocated in the Eulerian frame aligned fluid grid. Also, the manufactured variables can be defined at the final time to define error norms. Finally, the velocity components of interest are the ones defined in the Eulerian frame, and can be readily computed as
%
\begin{equation}
  \left\{ \begin{array}{c}  u \\ v \end{array} \right\}_p = [\mathbf{T}_{BE}]^T \left\{ \begin{array}{c}  u' \\ v' \end{array} \right\}_p
\end{equation}
%
We assume periodic boundary conditions on $\partial \Omega$. Then, the size of the domain $\Omega$ depends on the periodicity lengths for a rotated flow field.
When $\alpha=0$, a domain $\Omega$ of size $L_x=L_y=2\pi$ centered in $\mathbf{x}_{Ec}$ is sufficient for application of periodic boundary conditions. Then, for $\alpha=\arctan(1/r)$, $r=1,2,3,...$ the periodicity lengths on the Eulerian reference frame are $L_x=L_y=2\pi \sqrt{1+r^2}$.
That is, if $\alpha=\arctan(1)$,  the periodicity length increases to $L_x=L_y=\sqrt{8} \pi$. If $\alpha=\arctan(1/2)$, $L_x=L_y=\sqrt{20} \pi$.

% Figures showing each case with the geom and velocity field to show rotation:


% Figures rotcube_cc:
\begin{figure}[ht]
\centering
\includegraphics[trim = 5mm 2mm 7mm 2mm, clip, width=.365\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/rotated_cube_0deg_exp_mms_connvergence.pdf}
\includegraphics[trim = 30mm 2mm 7mm 2mm, clip, width=.305\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/rotated_cube_0deg_imp_mms_connvergence.pdf}
\includegraphics[trim = 30mm 2mm 7mm 2mm, clip, width=.305\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/rotated_cube_0deg_obs_mms_connvergence.pdf}
\caption[The {\ct Rotated Cube CC} accuracy order test case]{Zero degree rotation. The $L_2$ error defined on the whole domain at time $t = 0.95$ s is plotted as a function of grid spacing for $N=\{32,64,128,256,384\}$ points in each direction.}\label{fig:rotcube_cc_0deg_accuracy_order}
\end{figure}


\begin{figure}[ht]
\centering
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/rotated_cube_27deg_exp_mms_connvergence.pdf}
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/rotated_cube_27deg_imp_mms_connvergence.pdf}
\caption[The {\ct Rotated Cube CC} accuracy order test case]{$\simeq27$ degree rotation. The $L_2$ error defined on the whole domain at time $t = 0.95$ s is plotted as a function of grid spacing for $N=\{32,64,128,256,384\}$ points in each direction.}\label{fig:rotcube_cc_27deg_accuracy_order}
\end{figure}


Although the different orientation angles correspond to different periodicity lengths and domain sizes all cases were run with $N=\{32,64,128,256,384\}$ points in each direction for accuracy testing. The different inclination angles and grids were run with explicit and implicit time integration for scalars in the cut-cell region. The cases with zero inclination were also run with the $\&$OBST functionality for comparison. Calculations were run at constant $CFL=0.1$ and errors were computed on the solution at the first time level after $t=0.95$. The cases input files are found on \texttt{/Complex\_Geometry/rotated\_cube\_*deg\_*\_*.fds}. 

Overall the error norm behaves as $\mathcal{O}(\Delta x)$ for the scalar transport equation in cases where the cube is inclined. This is expected as the diffusive flux differences used are of this spatial accuracy. Note that for the non inclined cases of figure~\ref{fig:rotcube_cc_0deg_accuracy_order} the scalar $Z$ presents better accuracy, simply because the two point centroid to centroid difference coincides with a central difference. Note that the Godunov flux limited interpolation does not have a detrimental effect in accuracy, simply because advective fluxes reduce to zero at the boundary. 
The $\mathcal{O}(\Delta x^2)$ accuracy in $Z$ is also seen for $\&$OBSTs in the right figure of~\ref{fig:rotcube_cc_0deg_accuracy_order}, which use central differences for diffusive fluxes. 
For the $\&$OBSTs where the momentum equations are solved on velocity points next to the surface, velocities behave as $\mathcal{O}(\Delta x^2)$, in accordance to the discretization schemes used (right figure~\ref{fig:rotcube_cc_0deg_accuracy_order}). The immersed boundary velocity reconstruction on cut-faces provides spatial accuracies between $\mathcal{O}(\Delta x)$ and $\mathcal{O}(\Delta x^2)$ in $L_2$ norm across the different inclination angles, which are inline with what is expected for immersed boundary methods of this type. Note that pressure spatial accuracy is $\mathcal{O}(\Delta x)$ across all tests, consistent with the fractional step projection method used by FDS.    

% Figures for accuracy 27deg:
\begin{figure}[ht]
\centering
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/rotated_cube_45deg_exp_mms_connvergence.pdf}
\includegraphics[width=.49\textwidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/rotated_cube_45deg_imp_mms_connvergence.pdf}
\caption[The {\ct Rotated Cube CC} accuracy order test case]{$45$ degree rotation. The $L_2$ error defined on the whole domain at time $t = 0.95$ s is plotted as a function of grid spacing for $N=\{32,64,128,256,384\}$ points in each direction.}\label{fig:rotcube_cc_45deg_accuracy_order}
\end{figure}




\subsection{Poiseuille flow in channel created by two GEOMs}

Consider a gas with density $\rho=1.165$ kg/m$^3$, and molecular viscosity $\mu=0.025$ kg/m/s, being advected through a 2D channel of height $H=1$ m and length $L=10$ m, by means of a constant force field along the channel direction $F_x=1$ N/m$^3$. This force field is equivalent to a pressure drop $dp/dx=-1$ Pa/m. The analytical solution for the friction factor in this laminar flow is:
%
\begin{equation}
   f_{an} = 24/Re_H
\end{equation}
%
where in our case $Re_H= \rho \overline{u} H/ \mu \simeq 155$, being the analytical value for the bulk velocity $\overline{u}=3.\hat{3}$ m/s. The (Moody) friction factor formula for pressure losses is:
%
\begin{equation}
   \Delta p / L = f \frac{L}{H}  \frac{1}{2} \rho \overline{u}^2
\end{equation}
%
which can be used to compute the simulation error in $f$, as a function of the mean velocity evaluated on the numerical mesh. To evaluate the solution quality of the CC-IBM scheme we consider three sets of test cases. See figure~\ref{Fig:PoiseSketch}.
%
% figure
\begin{figure}[h]
      \centering
      \includegraphics[trim = 60mm 34mm 60mm 42mm, clip,
       width=0.75\linewidth]{../../../fig/fds/Geometry_Figures/Poiseuille_cc_Sketch.png}
      \caption{Poiseuille flow created by two geometries $\Omega_1, \Omega_2$: Location of channel walls respect to Cartesian grid lines (in blue).}
	\label{Fig:PoiseSketch}
\end{figure}
%
We define as $h$ the distance from the lower geometry wall to the closest grid line within it. The we consider discretizations with ${10,20,40,80}$ cells along the $z$ direction and cases: grid aligned $h=0$, grid dependent $h=\Delta z /3$, and fixed $h=\Delta z_{10}/11$, where $\Delta z_{10}$ is the wall normal cell size in the coarsest discretization used.

Friction factor errors computed as $|f-f_{an}|$ are shown in figures~\ref{Fig:PoiseConvg}a,b. Second order convergence of this error measure in all cases can be seen.
There is more noise on the fixed geometry case $h=\Delta z_{10}/11$, because the cut-cell size ratios and interpolation coefficients are changing from grid to grid. These error fluctuations are not seen when $h=\Delta z /3$, and cut-cells keep same sizes and interpolation coefficients respect to the Cartesian grids.
%
% figure
\begin{figure}[h]
      \centering
      \includegraphics[trim = 5mm 0mm 5mm 0mm, clip,
       width=0.49\linewidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/geom_poiseuille_convergence_theta0a.pdf}
      \includegraphics[trim = 5mm 0mm 5mm 0mm, clip,
       width=0.49\linewidth]{../FDS_Verification_Guide/SCRIPT_FIGURES/geom_poiseuille_convergence_theta0na.pdf}
      \caption{Poiseuille flow created by two geometries: Convergence trends on friction factor (Moody) for $h=0$, grid dependent $h=\Delta z/3$ and fixed $h=\Delta z_{10}/11$. Fluctuations on the friction factor error in the last case are produced by changes in cut-cell aspect ratios and velocity interpolation stencils as grids are refined.}
      \label{Fig:PoiseConvg}
\end{figure}
%


\subsection{Conservation test of isothermal helium plume around sphere}






%\subsection{Flow around Sphere at different Reynolds numbers}
%
%We look first at the canonical case of isothermal flow around a sphere of diameter $D$, for Reynolds numbers up to $Re=300$. This flow is steady for Reynolds numbers up to $Re \simeq 270$ where it transitions to an oscillatory state.
%
%The first set of steady cases are defined in a domain of size $[-3.7D,12.5D]$ in the streamwise $x$ direction, and $[-2.7D,2.7D]$ in the cross stream directions $y$ and $z$. A resolution of 40 cubic cells in $D$ is defined around the sphere and on the wake region up to $4.4D$. A coarser resolution of 20 cells in $D$ is defined in the far field.
%In figure~\ref{}a, the streamwise velocity curves along the axis of the sphere and from its center downstream are shown. The zero crossings of these curves on the fluid region define the length $x_s$ of the recirculation bubble behind the sphere.
%In figure~\ref{}b, $x_s$ is compared against results reported by Johnson and Patel~\cite{}. The agreement is very good, with a difference of the order of $10\%$ at $Re=50$ and $5\%$ on the other $Re=100,150,200$ data points.
%
%% Figs steady:
%%\begin{figure}
%%\includegraphics[width=3.2in]{SCRIPT_FIGURES/ }
%%\includegraphics[width=3.2in]{SCRIPT_FIGURES/ }
%%\caption{}
%%\end{figure}
%
%% Something on the force coefficients for steady.
%
%For the case of $Re=300$, the domain size defined is $[4D,14D]$ in the streamwise $x$ direction, and $[-3D,3D]$ in the cross stream directions. The blockage ratio for the sphere is in this case $2.18\%$. Tomboulides and Orzag~\cite{} provide an estimated laminar boundary layer thickness $\delta \simeq 0.065D$. We define a marginal resolution of $60$ cubic cells on $D$ on the sphere surface and in the wake region up to a downstream distance of $11D$. The Strouhal number from the streamwise velocity signal at $x=3.93D$ computed at oscillatory state is XXX, in excellent agreement with the the value of $St=0.136$ and $0.137$ provided in references~\cite{},\cite{} (computed from force coefficients). In figures~\ref{}a,b, mean and RMS streamwise velocity as function of axis positions downstream of the sphere are shown. Close agreement with reference~\cite{} is found. Lower RMS peaks on the near wake are noted on the current computation, likely consequence of a marginal resolution behind the sphere.
%
%
%% Figs Re=300:
%%\begin{figure}
%%\includegraphics[width=3.2in]{SCRIPT_FIGURES/ }
%%\includegraphics[width=3.2in]{SCRIPT_FIGURES/ }
%%\caption{}
%%\end{figure}
%
%% Something on the force coefficients for Re=300.


\subsection{Consistency test for OBST and GEOM of flow around heated cubes}


\subsection{Radiating Sphere}
\label{sphere_radiate}

A sphere with a radius of 1~m and a surface temperature of $T=773.15$~K is centered in a box that is 4~m on a side with a surface temperature of 0~K. The emissivity of all surfaces is 1, and there is no convective heat transfer to or from any surface. The heat flux from the hot sphere to the nearest cold wall surface 1~m away is expected to be $\sigma T^4/4=5.065$~kW/m$^2$. The results are shown in Fig.~\ref{sphere_radiate_fig}.
\begin{figure}[ht]
\includegraphics[width=3.2in]{SCRIPT_FIGURES/sphere_radiate.png}
\includegraphics[width=3.2in]{SCRIPT_FIGURES/sphere_radiate.pdf}
\caption[The {\ct sphere\_radiate} test case]{Left: a radiating hot sphere. Right: heat flux to the nearest wall surface.}
\label{sphere_radiate_fig}
\end{figure}


\subsection{Leaking Sphere}
\label{sphere_leak}

A hollow sphere with an inner radius of 1~m has a small cube within it that is blowing air into the sphere at a rate of $\dot{V}=0.008$~m$^3$/s. The sphere is assigned a leakage area of $A=0.0001$~m$^2$. It is expected that the pressure rise within the sphere shall be
\be
   \Delta p = \frac{\rho_\infty}{2} \left( \frac{\dot{V}}{A} \right)^2 \approx 3840 \; \hbox{Pa}
\ee   
as shown in Fig.~\ref{sphere_leak_fig}.
\begin{figure}[ht]
\centering
\includegraphics[width=3.2in]{SCRIPT_FIGURES/sphere_leak}
\caption[The {\ct sphere\_leak} test case]{Pressure rise in a leaking sphere.}
\label{sphere_leak_fig}
\end{figure}

%\subsection{Flow around heated curved pipe or something like that}


\clearpage

\section{Validation Tests}

\subsection{Vettori fire on sloped ceiling compartment experiments}





\subsection{Lattimer flow around inclined wall experiments}





\subsection{Askervein Hill}






\appendix

\section{Unstructured Geometry integral conservation tests}
\label{info:constests}

Consider an immersed body $B$, with surface defined by $N_\tau$ triangles, completely immersed in a domain $\Omega$ divided in $M$ cartesian mesh blocks. 
Three geometry conservation tests can be defined for the polyhedra definition around this body:
%
\begin{enumerate}
    
    \item Surface area: 
    The total body $B$ surface area can be computed from the triangulation $A_{B,\tau}=\sum_{i_\tau=1}^{N_\tau}{A_{i_\tau}}$, where $N_\tau$ is the number of triangles in the surface and $A_{i_\tau}$ their areas. It can also be calculated from the boundary cut-face areas computed by the cut-cell definition algorithm,  $A_{B,CF}=\sum_{icf=1}^{N_{bcf}}{A_{icf}}$, where $N_{bcf}$ is the number of cut-faces that have been computed on the body surface. These two areas should match.
    \newline
    
    \item Gas ($G$) region volume: 
    This is the volume of region $\Omega-\Omega_B$. It can be computed from the solid volume $V_B$ and the domain volume $V_\Omega$ as $V_{GB}=V_\Omega-V_B$. Here:
    %
    \begin{equation}
       V_\Omega=\sum_{nm=1}^M {V_{nm}} \: \: ; \: \: V_{nm}=L_{x,nm} L_{y,nm} L_{z,nm}
    \end{equation}
    %
    where $L_{x,nm}$, $L_{y,nm}$ and $L_{z,nm}$ are the $nm$ mesh sizes in each cartesian direction ($M=1$, $L_{x,1}=3D_B$, $L_{y,1}=3D_B$ and $L_{z,1}=8D_B$ in the cone problem of section~\ref{cone_test}). The discrete volume of the solid can be computed from the triangulation applying the divergence theorem to a field $\mathbf{F}_v=x \mathbf{\hat{i}}$, such that $\nabla \cdot \mathbf{\hat{F}}_v=1$, to the surface triangles:
    %
    \begin{equation}
      V_B = \int_{\Omega _{B}} \nabla \cdot \mathbf{F}_v \; dV = 
      \int_{\partial \Omega _{B}}  \mathbf{F}_v \cdot \hat{n} \; dA =
      \sum_{i=1}^{N_\tau}{(\mathbf{\hat{i}} \cdot \mathbf{\hat{n}}_i) x_{ci} A_i}
      \label{eq:Vb}
    \end{equation}
    %
    where $x_{ci}$ is the $i$th triangle centroid position in the $x$ direction, and $A_i$ its area.
    The $G$ region volume can also be computed from the sum of cut-cells and regular $G$ cells volumes over all $M$ cartesian meshes:
    %
    \begin{equation}
       V_{GRC} = \sum_{icc=1}^{N_{cc}}{V_{icc}} + N_{rg} \Delta x \Delta y \Delta z
       \label{}
    \end{equation}
    %
    where $N_{rg}$ the number of regular cells in the $G$ region, and $\Delta x,\Delta y,\Delta z$ the uniform grid cell sizes. The  $N_{cc}$ cut-cell volumes $V_{icc}$ are computed in similar manner to equation~\eqref{eq:Vb}:
    %
    \begin{equation}
    V_{icc} = \int_{\partial \Omega _{icc}}  \mathbf{F}_v \cdot \hat{n} \; dA = 
    \sum_{icf=1}^{ncf} (\mathbf{\hat{i}} \cdot \mathbf{\hat{n}}_{icf}) \int_{\Omega _{icf}} x \; dA
    \label{eq:vol2}
    \end{equation}
    %
    where $ncf$ is the number of faces (cut-face or regular) which are boundary of cut-cell $icc$.
    
    In this test $V_{GB}$ and $V_{GRC}$ should agree.
    \newline 
    
    \item $G$ region centroid location:
    The centroid of the $G$ region can be computed from cut-cells and regular $G$ cells as:
    %
    \begin{equation}
       \mathbf{x}_{GRC} = \frac{1}{V_{GRC}} 
       \left[  \sum_{icc=1}^{N_{cc}}{\mathbf{x}_{c(icc)}V_{icc}} +  \sum_{irg=1}^{N_{rg}}{\mathbf{x}_{c(irg)} \Delta x \Delta y \Delta z} \right] 
    \end{equation}
    %
    Where $\mathbf{x}_{c(irg)}$ is the cell center location for regular cell $irg$, and the scalar components in each coordinate direction  $\{\mathbf{\hat{i}},\mathbf{\hat{j}},\mathbf{\hat{k}}\}$ of cut-cell centroids $\mathbf{x}_{c(icc)}$ are obtained from:
    %
    \begin{equation}
    x_{ci(icc)} = \frac{1}{V_{icc}} \int_{\Omega _{icc}} \nabla \cdot \mathbf{F}_{c(i)} \; dV =
    \frac{1}{V_{icc}} \int_{\partial \Omega _{icc}} \mathbf{F}_{c(i)} \cdot \hat{n} \; dA = 
  \frac{1}{2V_{icc}} \sum_{icf=1}^{ncf} (\mathbf{\hat{e}}_i \cdot \mathbf{\hat{n}}_{icf}) \int_{\Omega _{icf}} x_i^2 \; dA 
    \label{eq:xcen}
    \end{equation}
    %
    where $\mathbf{F}_{c(i)}=1/2 x_i^2 \mathbf{\hat{e}}_i$. Here, $x_i=\{x,y,z\}$ and $\mathbf{\hat{e}}_i=\{\mathbf{\hat{i}},\mathbf{\hat{j}},\mathbf{\hat{k}}\}$.
    
    The centroid of the domain $\Omega$ is:
    \begin{equation}
        \mathbf{x}_\Omega = \frac{1}{V_\Omega} \sum_{nm=1}^M  \mathbf{x}_{nm} V_{nm}
    \end{equation}
    where $\mathbf{x}_{nm}=\mathbf{x}_{low,nm}+\frac{1}{2} \left( L_{x,nm},L_{y,nm},L_{z,nm} \right)$, where $\mathbf{x}_{low,nm}$ is the vector of low boundary coordinates for mesh $nm$. In the single mesh cone problem (section~\ref{cone_test}) $\mathbf{x}_{low,1}= \left( -1.5D_B,-1.5D_B,-2D_B \right)$.
    
    Also, the centroid $\mathbf{x}_B$ of the discretized solid can be computed applying equation~\eqref{eq:xcen} to the surface triangulation $\tau$, in same manner to what was done for cut-faces. From domain and solid the centroid of the $G$ region can be obtained as:
    %
    \begin{equation}
      \mathbf{x}_{GB}=\frac{1}{V_{GB}} \left[ \mathbf{x}_\Omega V_\Omega - \mathbf{x}_B V_B \right]
    \end{equation}
    %
    Consistency of the computed centroids requires that $\mathbf{x}_{GRC}$ and $\mathbf{x}_{GB}$ should correspond.
    
\end{enumerate}



%\paragraph{Abstract} In this chapter, we discuss implementation of a hybrid cutcell direct-forcing immersed boundary method (CIBM) to handle complex geometry in the Fire Dynamics Simulator (FDS).  Particular attention is paid to the challenges associated with thin (zero thickness) obstructions in variable density flows with heat and mass transfer.  The triangular facet is chosen as the basic element for construction of complex geometry.  Reasons behind this choice are discussed.  A linked-list data structure for the facets is presented.  One-sided differencing and methods for overcoming numerical diffusion issues are discussed.  Verification cases for momentum, heat, and mass transport are presented.
%
%\subsection{Introduction}
%
%The Fire Dynamics Simulator (FDS) has always used a simple direct forcing immersed boundary method (IBM) to handle flow obstructions.  Generally, obstructions are restricted to conform to the underlying Cartesian mesh, leading to stair-step approximations to curvilinear geometries (see Fig.~\ref{fig:circle}).
%
%\begin{figure}
%\begin{center}
%\begin{tabular}{cc}
%\begin{tikzpicture}
%\draw[line width=1mm,fill=lightgray] (1.8,1.8) circle [radius=1.5];
%\draw[help lines] (0,0) grid (5,5);
%\draw[ultra thick,->] (2.66,4.5) -- (3.33,4.5);
%\draw[ultra thick,->] (2.66,3.5) -- (3.33,3.5);
%\draw[ultra thick,->] (3.66,4.5) -- (4.33,4.5);
%\draw[ultra thick,->] (3.66,3.5) -- (4.33,3.5);
%\draw[thick] (2.665,3.0255) -- (3.335,3.9745);
%\draw[fill] (3.335,3.9745) circle [radius=0.1];
%\end{tikzpicture}
%\hspace{1cm}&\hspace{1cm}
%\begin{tikzpicture}
%\draw[line width=1mm,fill=lightgray] (1,0) to (1,1) to (0,1) to (0,2) to (1,2) to (1,3) to (3,3) to (3,1) to (2,1) to (2,0) to (.95,0);
%\draw[help lines] (0,0) grid (5,5);
%\end{tikzpicture}
%\end{tabular}
%\caption{(Left) Desired geometry (circle). (Right) FDS representation---a cell is solid if more than half its volume is occupied by the desired geometry.}
%\label{fig:circle}
%\end{center}
%\end{figure}
%
%
%There are two shortcomings of the block geometry method.  First, it is zeroth-order accurate for curvilinear geometries \cite{Fadlun:temp}.
%In other words, as the grid is refined---even after laborious reconstruction of the block geometry---we never recover the exact solution.  Second, creation of even simple geometries can be tedious. For example, users will go to incredible lengths to carve out rounded tunnels for circular ducts, a task that is fairly simple for most unstructured codes (not that the meshes tend to be much better!).
%
%In direct forcing IBM, the force term in the momentum equation is set to drive a specific velocity component to a desired value.  In the block geometry approach (Fig.~\ref{fig:circle}, right), the impermeability condition is enforced by setting the normal component of velocity to zero on the solid surface. But as shown by Fadlun \cite{Fadlun:temp}, it is relatively simple to obtain a second-order solution for the velocity field: the velocity component is found from a linear interpolation between the velocity of the solid surface (no slip, impermeability) and a value at some distance normal to the surface.  This is illustrated on the left in Fig.~\ref{fig:circle}.  Away from the surface the staggered $u$ velocity components are interpolated to the point marked by the solid dot, which is twice the normal distance from the surface as the nearest velocity component.
%
%More elaborate schemes exist to obtain the velocity near the surface \cite{Balaras:temp,Choi:temp,McDermott:temp}.  But the emphasis here is on the implementation of a scheme to transport heat and mass near thin obstructions (an issue yet to be addressed in the literature).  Therefore, Fadlun's simple interpolation method is sufficient for momentum.   In 2D, complex surfaces are constructed from a series of connected line segments. In 3D, surfaces are constructed from triangular facets.  Our goal is to minimize numerical diffusion across thin obstructions arbitrarily oriented relative to the Cartesian mesh.  The cause of this diffusion is discussed in more detail below.
%
%In the following sections, we first discuss the linked-list data structure used to connect the surface elements to the grid, similar to the wall cell type used in the basic FDS scheme.  Next, we illustrate the momentum scheme for both simple and complex flow geometries.  We then discuss the unique issues faced by the FDS flow solver for simulating variable density flows with thin obstructions not aligned with the mesh.  Hopefully, we can develop a scheme that handles these problems!
%
%\subsection{Data Structure}
%
%\subsubsection{Block Geometry}
%
%In the basic FDS method, a wall cell (WC) is one face of a gas phase cell.  A WC derived type stores the indices of the gas phase cell\footnote{Here I am not getting into the difference between wall cell and 1D in-depth derived types.  When I say WC, in most cases I mean precisely {\ct WC\%ONE\_D}.}.  It also stores face values of density, species mass fraction, and temperature, and the wall-normal cell dimension used to construct gradients.  The relevant variables in the code are:
%
%\begin{table}[h!]
%\begin{tabular}{ll}
%{\ct WC\%IIG, WC\%JJG, WC\%KKG} & gas phase cell indices \\
%{\ct WC\%UW}                    & wall-normal velocity component \\
%{\ct WC\%RHO\_F}                & gas phase mass density at the wall \\
%{\ct WC\%ZZ\_F}                 & lumped species mass fraction at the wall \\
%{\ct WC\%TMP\_F}                & temperature at the wall \\
%{\ct WC\%DN}                    & wall-normal cell dimension \\
%{\ct WC\%IOR}                   & orientation from wall to cell center
%\end{tabular}
%\end{table}
%
%The face values serve as Dirichlet boundary conditions.  If Neumann boundary conditions are needed, they are constructed from a one-sided difference.  For example, the temperature gradient at the wall is obtained from
%\begin{equation}
%\frac{\partial T}{\partial n} \approx \frac{T_g - T_f}{\delta n/2}
%\end{equation}
%
%\subsubsection{Complex Geometry}


%As we move toward complex shapes immersed within the Cartesian mesh, to the extent possible we have tried to maintain a correspondence with the WC data type.  The analog in the new approach is, in 3D, a triangular facet, FC. This choice was made for the following reasons:
%\begin{enumerate}[{\,\,\,\,(}a{)}]
%\item For simplicity, we wanted to restrict ourselves to one generic element type.
%\item Curvilinear surfaces can be approximated to second-order accuracy with linear elements.
%\item Three points are unambiguously planar.
%\item Algorithms for triangle/box intersection are readily available.
%\end{enumerate}
%
%Note that a basic WC face must be constructed from two FC faces (see Fig.~\ref{fig:facetcell}), thus doubling computational effort for the same resolution.  The FC approach makes up for this deficiency with the following benefits:
%\begin{enumerate}[{(}i{)}]
%\item We have the flexibility to create arbitrarily complex shapes.
%\item The resolution of the solid surface need not be the same as the gas phase.  When coarser resolution is permissible, we achieve a cost savings.
%\end{enumerate}
%
%%The first item may be obvious, but the second deserves elaboration. For many objects in the scenario, the boundary conditions may not change much or may not be important enough to finely resolve.  Suppose there is a door that is only partially open (say 45$^\circ$ angle) and the gas phase resolution within the room is 10 cm.  A standard doorway in the U.S. is 91 cm by 203 cm.  Therefore, resolving each gas phase cell with a wall cell on the door would require roughly $9 \times 20 \times 2 = 360$ facets.  But if detailed resolution of the door is unnecessary, we can get away with only 2 facets (see Fig.~\ref{fig:opendoor}).
%
%\begin{figure}
%\begin{center}
%\begin{tikzpicture}
%
%\draw (0,0)--(0,3);
%\draw (0,3)--(3,3);
%\draw (3,3)--(3,0);
%\draw (3,0)--(0,0);
%
%\draw (1,4)--(4,4);
%\draw (4,4)--(4,1);
%
%\draw (0,3)--(1,4);
%\draw (3,3)--(4,4);
%\draw (3,0)--(4,1);
%
%\draw (0,0)--(3,3);
%\draw (3,3)--(4,1);
%\draw (1,4)--(3,3);
%
%\end{tikzpicture}
%
%\caption{Triangular facets covering a cubic cell.  Two facets are needed for each face, as opposed to one wall cell.}
%\label{fig:facetcell}
%\end{center}
%\end{figure}
%
%
%%\begin{figure}
%%\begin{center}
%%\begin{tikzpicture}[scale=0.75]
%%
%%\draw (0,6) -- (2,8);
%%\draw (0,0) -- (2,2);
%%
%%\draw (12,6) -- (10,8);
%%\draw (12,0) -- (10,2);
%%
%%\draw (2,8) -- (10,8);
%%\draw (2,8) -- (2,2);
%%\draw (10,8) -- (10,2);
%%
%%\draw (4,6.5) -- (6.5,6.5);
%%\draw (4,2) -- (4,6.5);
%%\draw (2,2) -- (4,2);
%%\draw (6.5,2) -- (10,2);
%%
%%\draw (0,0) -- (0,6);
%%\draw (12,0) -- (12,6);
%%\draw (0,0) -- (12,0);
%%
%%\draw[help lines] (1.499,1.499) grid [step=.25cm] (10.5,7.5);
%%
%%\draw[help lines,fill=white] (5,1) to (5,5.5) to (6.5,2) to (5,1);
%%\draw[help lines,fill=white] (6.5,6.5) to (5,5.5) to (6.5,2) to (6.5,6.5);
%%
%%\draw (6.5,2) -- (6.5,6.5);
%%\draw (6.5,6.5) -- (5,5.5);
%%\draw (6.5,2) -- (5,1);
%%\draw (5,1) -- (5,5.5);
%%\draw (5,5.5) -- (6.5,2);
%%
%%\draw[help lines] (5.7499,1.5) grid [step=.25cm] (7.5,7.5);
%%
%%\end{tikzpicture}
%%
%%\caption{Partially open door composed of two facets, illustrating that the solid resolution may be coarser than the gas phase resolution.}
%%\label{fig:opendoor}
%%\end{center}
%%\end{figure}
%
%
%\subsection{Momentum}
%\label{sec:momentum}
%
%Consider the following semi-implicit (fractional step) update for the $i$th component of momentum:
%\begin{equation}
%\label{eqn:mom}
%\frac{u_i^{n+1} - u_i^n}{\delta t} = -\left(F_i^n + \frac{\delta H^{n+1}}{\delta x_i}\right)
%\end{equation}
%This is a fractional step method because a Poisson equation is solved to obtain the pressure term, $H$ (not discussed here). In a direct-forcing immersed boundary method, the force term is replaced by
%\begin{equation}
%\label{eqn:force}
%F_i^* = - \left( \frac{u_i^{ibm} - u_i^n}{\delta t} + \frac{\delta H^*}{\delta x_i}\right)
%\end{equation}
%which drives the velocity component $u_i^{n+1}$ toward $u_i^{ibm}$.  The superscript asterisk indicates that the initial guess for the pressure term is taken from the previous time step.  The update may be iterated to drive the error between $u_i^{n+1}$ and $u_i^{ibm}$ to a desired tolerance.  Note, however, that when this approach is embedded in a second-order predictor-corrector scheme the error in most cases is remarkably small when $H^n$ is used in the force term (no iteration).
%
%\subsubsection{Block Geometry}
%
%For the basic block geometry method used in FDS, the impermeability condition at a surface is enforced by setting $u_i^{ibm} = 0$. If mass transfer is present the surface-normal IBM velocity component is set to $u_i^{ibm} = u_n$, where $u_n$ is the wall-normal component of velocity satisfying the scalar transport equation.  Note that $u_n$ is stored as {\ct WC\%UW} in the code.
%
%\subsubsection{Complex Geometry}
%
%For complex geometry the surface usually does not align itself with the face of a Cartesian grid cell. A given velocity component, therefore, is not on a boundary, but must still satisfy the governing equations subject to the prescribed boundary conditions (position and properties of the surface).  The essential task of a direct-forcing immersed boundary method is to set $u_i^{ibm}$ to achieve a stable and accurate numerical solution as well as an adequate representation of the local flow physics (surface stress, pressure gradient, etc.) without the luxury of infinitely high resolution.
%
%In FDS, there are three IBM options for momentum, set by, for example,
%\small
%\begin{verbatim}
%&MISC IMMERSED_BOUNDARY_METHOD=1/ ! options are 0,1,2
%\end{verbatim}
%\normalsize
%Following is a brief description of each option. The basic linear interpolation method of Fadlun \cite{Fadlun:temp} is given by option (1), and will be used for the examples to follow later in this paper. Refer to the diagram in Fig.~\ref{fig:position_definitions} for position definitions.
%
%\paragraph{\ct IMMERSED\_BOUNDARY\_METHOD=0}
%
%The velocity component force term is unaltered if it falls on the positive side of the surface.  Let $\mathbf{x}_{velo}$ denote the location of the velocity component, $\mathbf{v}_i, (i=1,2,3)$ are the vertices of a face, $\mathbf{n} = (\mathbf{v}_2-\mathbf{v}_1) \times (\mathbf{v}_3-\mathbf{v}_1)$ is the face normal, and let $\mathbf{r} \equiv \mathbf{x}_{velo} - \mathbf{v}_1$ (any of the vertices will do); then if $\delta n =\mathbf{r} \cdot \mathbf{n} > 0$ the velocity component is processed as usual; if the dot product is negative the velocity is ``inside'' the obstruction and is set to zero.  This effectively leads to a zeroth-order (Lego block) representation of the surface.
%
%\paragraph{\ct IMMERSED\_BOUNDARY\_METHOD=1}
%
%The basic linear interpolation scheme of Fadlun starts by identifying a point $\mathbf{x}_{int}$ (see Fig.~\ref{fig:position_definitions}) which is twice the normal distance from the surface as the position of the velocity component, $\mathbf{x}_{velo}$.  The value $u_i(\mathbf{x}_{int}) = u_{i,int}$ is found from a trilinear interpolation of the neighboring staggered velocity component values from the previous time step.  The value of $u_i^{ibm}$ is then taken from a linear interpolation between $u_{i,int}$ and the wall value of the velocity component, $u_{i,w}$ (usually zero).  Thus, $u_i^{ibm} = \frac{1}{2}(u_{i,int} + u_{i,w})$.
%
%\paragraph{\ct IMMERSED\_BOUNDARY\_METHOD=2}
%
%Here the linear profile is replaced by a prescribed wall-normal profile for the mean streamwise velocity.  This requires that we first define a new streamwise coordinate system.  Referring to Fig.~\ref{fig:position_definitions}, let $\mathbf{e}_i, (i=1,2,3)$ define an orthonormal set of basis vectors for our Cartesian grid system and let $\mathbf{\bar{e}}_i, (i=1,2,3)$ define the new streamwise system where $\mathbf{\bar{e}}_1$ is the streamwise direction and $\mathbf{\bar{e}}_3$ is the wall-normal direction.  It is useful to define the new system so that the velocity component in $\mathbf{\bar{e}}_2$ is always zero.  This is achieved as follows: Define a tangent vector $\mathbf{t} \equiv \mathbf{n} \times (\mathbf{u} - \mathbf{u}_w)$ where $\mathbf{u}_w$ is the velocity of the wall (usually zero) and $\mathbf{n}$ is a wall-normal vector. The streamwise basis vector is then found from $\mathbf{\hat{s}} = \mathbf{\hat{t}} \times \mathbf{\hat{n}}$.  The hat represents a unit normal, e.g., $\mathbf{\hat{n}} \equiv \mathbf{n}/|\mathbf{n}|$.  The new streamwise system is then given by $\mathbf{\bar{e}_1}=\mathbf{\hat{s}}$, $\mathbf{\bar{e}_2}=\mathbf{\hat{t}}$, and $\mathbf{\bar{e}_3}=\mathbf{\hat{n}}$. (To conform to the convention used in the atmospheric community, in FDS, the ``vertical'' or surface-normal component is taken as the ``$z$'' component.  Hence, 2D flows are in the $x-z$ plane.)  The elements of the directional cosine matrix are now simply $a_{ij} = \mathbf{e}_i \cdot \mathbf{\bar{e}}_j$.  Velocity vectors may be transformed from the grid system to the streamwise system via $\bar{u}_k = a_{ik} u_i$ and back again via $u_i = a_{ik} \bar{u}_k$ (see, e.g., \cite{Pope:2000} Appendix B).
%
%To evaluate the streamwise profile, all three staggered velocity components are interpolated to $\mathbf{x}_{int}$ and the vector is transformed into the streamwise system.  The streamwise component at $\mathbf{x}_{int}$ is used to define the profile $f(n)$, which is then used to specify the streamwise component at the position of interest ($\mathbf{x}_{velo}$ in the grid system).  Thus, we obtain $\bar{u}_1^{ibm} = f(\delta n)$.  Here note that the bar indicates we are in the streamwise coordinate system.  We need to specify the other two components before we can transform back to the grid system.  As mentioned, the tangential component is zero by construction. Specification of the normal component is currently an active area of research \cite{Choi:2007:temp}.  What complicates matters is that we have two invariants to satisfy: (1) the divergence and (2) the kinetic energy. The kinetic energy is the easier of the two because it does not involve differences.  Still, since all we know from the model is a mean profile, some sort of \emph{ad hoc} assumption must be made.  Given $k_{int}=\frac{1}{2}u_i u_i$ at $\mathbf{x}_{int}$, we use a linear interpolation to determine $k(\delta n) = \frac{1}{2}(k_{int} + k_w)$, where $k_w$ is the kinetic energy at the wall (usually zero).  The magnitude of the normal component is then found from $|\bar{u}_3^{ibm}| = +\sqrt{2k(\delta n) - (\bar{u}_1^{ibm})^2}$.  Finally, the component value of interest is obtained from the transformation $u_i^{ibm} = a_{ik} \bar{u}_k^{ibm}$, where only the $i$th component is retained at the staggered location.
%
%\begin{figure}
%\begin{center}
%\begin{tikzpicture}[scale=1.25]
%
%\draw [fill=lightgray] (0.2,0.5) -- (2.2,3.8) -- (4.5,3.1) -- (0.2,0.5);
%
%\draw[help lines] (-3,-1) grid [step=2cm] (5,7);
%
%\draw (0.30,0.25) node {$\mathbf{v}_1$};
%\draw (4.80,3.10) node {$\mathbf{v}_2$};
%\draw (2.20,4.00) node {$\mathbf{v}_3$};
%
%\draw [->, very thick, >=latex] (-2.5,5.0) -- (-1.5,5.0);
%\draw [->, very thick, >=latex] (-2.5,3.0) -- (-1.5,3.0);
%\draw [->, very thick, >=latex] (-0.5,5.0) -- (0.5,5.0);
%\draw [->, very thick, >=latex] (-0.5,3.0) -- (0.5,3.0);
%\draw (0.5,3.3) node {$\mathbf{x}_{velo}$};
%
%\draw [->, very thick, >=latex] (0.2,0.5) -- (0.0,3.0);
%\draw (-0.15,1.8) node {$\mathbf{r}$};
%
%\draw [->, very thick, >=latex] (0.2,0.5) -- (-1,2.0);
%\draw (-0.5,1.0) node {$\mathbf{n}$};
%
%\draw [->, >=*] (1.1,1.545) -- (-1.3,4.7);
%\draw (-1.5,4.3) node {$\mathbf{x}_{int}$};
%
%\draw [|-] (0.9,1.37) -- (0.57,1.8);
%\draw [-|] (0.24,2.25) -- (-0.19,2.85);
%\draw (.4,2.02) node {$\delta n$};
%
%\draw [->, very thick, >=latex] (1.08,1.55) -- (2.35,1.55);
%\draw [->, very thick, >=latex] (1.1,1.545) -- (0.45,2.4);
%\draw [->, very thick, >=latex] (1.1,1.545) -- (1.75,2.6);
%\draw (2.9,1.55) node {$\mathbf{\hat{s}} \, (\mathbf{\bar{e}}_1)$};
%\draw (1.9,2.60) node {$\mathbf{\hat{t}}$};
%\draw (0.6,2.60) node {$\mathbf{\hat{n}}$};
%
%\draw [->, very thick, >=latex] (-2,0) -- (-2,1);
%\draw [->, very thick, >=latex] (-2,0) -- (-1,0);
%\draw (-0.9,-0.25) node {$\mathbf{e}_1$};
%\draw (-2.35,1) node {$\mathbf{e}_3$};
%
%\draw (1,4.7) .. controls (2.5,2) .. (1.1,1.545);
%\draw [->] (-1.15,4.500) -- (1.10,4.500);
%\draw [->] (-0.58,3.750) -- (1.52,3.750);
%\draw [->] ( 0.20,3.003) -- (1.90,3.003);
%\draw [->] ( 0.54,2.280) -- (2.15,2.280);
%\draw [->] ( 0.90,1.800) -- (1.80,1.800);
%\draw (1.5,5.0) node {$f(n)$};
%
%\end{tikzpicture}
%
%\caption{Position definitions for a single facet.  The function $f(n)$ represents the mean streamwise velocity profile in the wall-normal direction.}
%\label{fig:position_definitions}
%\end{center}
%\end{figure}
%
%\paragraph{Remark} As can be seen, in going from IBM methods 0 to 2 the level of complexity increases significantly.  For the remainder of this paper, we utilize method 1.  The main reason is that our focus here is on the additional complexities that arise with considering heat and mass transfer near infinitely thin obstructions.
%
%%\subsubsection{Convergence of IBM 1}
%%\label{sec:momentum_convergence}
%%
%%In this section we show second-order convergence for IBM 1 by examining the 2D flow near a rotating disk.
%
%\subsection{Heat and Mass Transfer}
%
%This is where the formulation gets tricky.  Both heat and mass transfer present challenges for a typical IBM formulation.  For heat transfer, numerical diffusion from the gas phase domain into the solid domain implies that global energy is not conserved.  For mass transfer, when considering the ``mixed is burnt'' model for combustion chemistry, any numerical diffusion of reactant species across a solid boundary can be fatal to a calculation.  For these reasons, we have chosen to implement a hybrid cutcell immersed boundary method (CIBM), where the cutcells on the immersed boundary surface are treated by an unstructured finite volume solver for heat and mass transport.  This ensures strict conservation of mass and energy in the fluid domain.
%
%\subsubsection{To-Do List}
%
%\begin{enumerate}
%\item Identify cutcells.  Find the Cartesian cells that intersect with solid tet volumes.  To handle zero thickness obstructions we either need the ability to treat facets by themselves or we may allow ``volumes'' with zero volume.  To be discussed.
%\item Compute cutcell (CC) volumes.  Let $V_{cell}$ denote the volume of the Cartesian cell.  Let $V_{tet}$ denote a tetrahedral volume.  The cutcell volume is $V_{cc} = V_{cell} - ( V_{cell} \cap V_{tet} )$.
%\item Merge small cells.  If $V_{cc} < V_{tol}$, merge with a neighboring cell (Cartesian cell or cutcell).  This limits CFL time step restrictions.
%\item Compute surface areas and unit normal vectors for final cutcell volumes.
%\item Interpolate scalars (density, mass fractions, temperature) to CC volume faces (see below).
%\item Interpolate velocity components to CC faces (see below).
%\item Update finite volumes equations (see below).
%\item Construct Cartesian divergence from CC divergence integral (see below).
%\item Update Poisson equation as usual (see FDS Tech Guide).
%\item Reconstruct velocity field to satisfy projection and CC constraints (see below).
%\end{enumerate}
%
%\newpage
%
%\subsubsection{FDS Grid Data Structure}
%
%
%In the following, variables and proposed data lists for cut-cells and faces to be used in the FDS CCIBM code are described. Data structures from the prototype Matlab code and linked list definitions from the current FDS IBM implementation are included for clarity. Names can be changed if there are issues with the existing FDS namespace, etc. Also, we keep it general for now.
%
%
%\subsubsection*{Some type defining and indexing parameters:}
%
%The following parameters are used in defining the type (or status) of cell or cell face on Cartesian grid or solid surfaces: \\
%\texttt{INTEGER, PARAMETER :: GASPHASE     = -1} \\
%\texttt{INTEGER, PARAMETER :: CUTCF         =   0} \\
%\texttt{INTEGER, PARAMETER :: SOLID              =   1} \\
%\texttt{INTEGER, PARAMETER :: INBOUNDARY =   2} \\
%
%
%
%\subsubsection*{New Cartesian grid variables:}
%
%- \textit{Cell-centered} variables for cell-type identification, unknown numbering, etc.:  \\
% \texttt{INTEGER} \texttt{CCVAR(I,J,K,ICVAR)}: Here,  \texttt{ICVAR} can be equal to:   \\
%%
%\begin{itemize}
%
%  \item \texttt{INTEGER, PARAMETER :: IGSC = 1} such that: \\
%           - \texttt{CCVAR(I,J,K,IGSC) == GASPHASE} is a regular \texttt{GASPHASE} cell. \\
%           - \texttt{CCVAR(I,J,K,IGSC) == CUTCF} is a cut-cell, i.e. cell transversed by the solid wet surface. \\
%           - \texttt{CCVAR(I,J,K,IGSC) == SOLID} is a regular cell immersed completely in the solid region.
%
%  \item \texttt{ICVAR} can take other values (2, 3,... to be added later) to map unknowns to \texttt{CUT\_CELL} data, linear systems, etc.
%
%
%\end{itemize}
%%
%- \textit{Face-centered} variables for face-type identification, generating data necessary for matrix and vector building (face based),  etc.:  \\
%\texttt{INTEGER} \texttt{FCVAR(I,J,K,IFVAR)}: Here,  \texttt{IFVAR} can be equal to:   \\
%%
%\begin{itemize}
%
%  \item \texttt{INTEGER, PARAMETER :: FGSC = 1} such that: \\
%           - \texttt{FCVAR(I,J,K,FGSC) == GASPHASE} is a regular \texttt{GASPHASE} face. \\
%           - \texttt{FCVAR(I,J,K,FGSC) == CUTCF} is a cut-face, i.e. face transversed by the solid wet surface. \\
%           - \texttt{FCVAR(I,J,K,FGSC) == SOLID} is a regular face immersed completely in the solid region.
%
%  \item \texttt{FCVAR} can take other values (2,3,... to be added later) to map face indexes to \texttt{CUT\_FACE} data, etc.
%
%
%\end{itemize}
%%
%These variables are defined by mesh, so they can be included as \texttt{ALLOCATABLE} in the \texttt{TYPE\_MESH}, and be allocated if using CCIBM.
%
%
%\subsubsection*{Data structures currently used on the Matlab prototype code:}
%
%- Data structure \texttt{CUT\_FACE}: Contains all \texttt{GASPHASE} and \texttt{INBOUNDARY} cut-faces. The latter are defined in the immersed solid surface and are used to impose boundary conditions for scalars. The total number of cut-faces is \texttt{NCUTFACE}. For a given cut-face \textit{icf}, the data structure has fields:
%%
%\begin{itemize}
%
%  \item  \texttt{CUT\_FACE(icf).XYZVERT}: $x,y,z$ coordinates of cut-face \textit{icf} vertex points. They are composed of intersection points between solid and Cartesian cell faces, and cell or facet vertices. These can be stored on a Global array as discussed with Glenn, and \texttt{XYZVERT} can be replaced by an index vector. They are used to define geometric data and also for plotting purposes (would be used in SMV).
%
%  \item  \texttt{CUT\_FACE(icf).XYZCEN}: $x,y,z$ coordinates of cut-face \textit{icf} area centroid.
%
%  \item  \texttt{CUT\_FACE(icf).STATUS}: can be \texttt{GASPHASE} or \texttt{INBOUNDARY}.
%
%  \item  \texttt{CUT\_FACE(icf).IJK}: If \texttt{CUT\_FACE(icf).STATUS == GASPHASE}, \texttt{ijk} denotes the Cartesian grid \textit{face} the cut-face belongs to. If  \texttt{CUT\_FACE(icf).STATUS == INBOUNDARY},  \texttt{ijk} denotes the Cartesian grid \textit{cell} this cut-face belongs to.
%
%  \item  \texttt{CUT\_FACE(icf).AXIS}: If \texttt{CUT\_FACE(icf).STATUS == GASPHASE}, this face is normal to which axis? \texttt{AXIS} can be \texttt{== IAXIS=1} for x-faces, \texttt{== JAXIS=2} for y-faces and \texttt{== KAXIS=3} for z-faces.  \texttt{CUT\_FACE(icf).STATUS == INBOUNDARY}, then  \texttt{AXIS} is zero.
%
%  \item \texttt{CUT\_FACE(icf).NIJK}: normal unit vector to the cut-face. If \texttt{CUT\_FACE(icf).AXIS == IAXIS} then \texttt{NIJK=[1. 0. 0.]}, if
%  \texttt{CUT\_FACE(icf).AXIS == JAXIS} then \texttt{NIJK=[0. 1. 0.]}, \texttt{CUT\_FACE(icf).AXIS == KAXIS} then  \texttt{NIJK=[0. 0. 1.]}. If \texttt{CUT\_FACE(icf).AXIS == 0},  \texttt{NIJK} is the normal outside of the body for that \texttt{INBOUNDARY} cut-face.
%
%  \item \texttt{CUT\_FACE(icf).BCTYPE(1:NSCALARS)}: boundary condition type, integer that for \texttt{INBOUNDARY} cut-faces, specifies the type of boundary condition imposed for scalars (i.e. constant mass flux, wall, etc.).
%
%  \item \texttt{CUT\_FACE(icf).BCVAL(1:NSCALARS)}: mean boundary condition value assigned for this \texttt{INBOUNDARY} cut-face, for each scalar.
%
%\end{itemize}
%%
%Some of the data here is redundant and won't need to be filled up. Now, as in this construct cut faces can be either \texttt{GASPHASE} or \texttt{INBOUNDARY}, we might want to separate these in two structs (i.e. CUT\_FACE\_GP and CUT\_FACE\_BD) to avoid conditionals on cut-face loops.
%
%
%- Data structure \texttt{CUT\_CELL}: Contains all \texttt{GASPHASE} cut-cells. The total number of cut-cells is \texttt{NCUTCELL}. For a given cut-cell \textit{icc}, the data structure has fields:
%%
%\begin{itemize}
%
%  \item \texttt{CUT\_CELL(icc).vol}: cut-cell volume.
%
%  \item \texttt{CUT\_CELL(icc).XYZCEN}:  $x,y,z$ coordinates of cut-cell \textit{icc} volume centroid.
%
%  \item \texttt{CUT\_CELL(icc).NCFACES}: number of faces in the \texttt{CUT\_FACE} data structure.
%
%  \item \texttt{CUT\_CELL(icc).CFACES}: index vector to location of cut-faces in \texttt{CUT\_FACE} data structure.
%
%  \item \texttt{CUT\_CELL(icc).IJK}: denotes the Cartesian grid \textit{cell} this cut-cell belongs to.
%
%  \item \texttt{CUT\_CELL(icc).HINTFLG}: Integer flag that denotes if the pressure $H$ needs to be interpolated to the \texttt{CUT\_CELL(icc).XYZCEN} from the \texttt{GASPHASE} surroundings, and type of interpolation to perform.
%
%\end{itemize}
%%
%Other fields will need to be incorporated to this data structure. Among them an unknown number for scalar and possibly the pressure H, if pressure is dealt with using the unstructured scheme.
%
%
%
%\subsubsection*{Data structures currently used on the FDS IBM code:}
%
%\begin{myfont}
%
%\noindent TYPE CUTCELL\_LINKED\_LIST\_TYPE \\
%\indent  INTEGER :: INDEX                                   ! data \\
%\indent  TYPE(CUTCELL\_LINKED\_LIST\_TYPE), POINTER :: NEXT    ! nxt el \\
%\indent   REAL(EB) :: AREA                                   ! cutcell area for index \\
%\noindent END TYPE CUTCELL\_LINKED\_LIST\_TYPE \\
%
%\noindent TYPE CUTCELL\_TYPE \\
%\indent   REAL(EB) :: VOL,RHO,TMP,DIV,A(6),S,N(3) \\
%\indent   REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ \\
%\indent   ! identify each face with area contribution to current cutcell \\
%\indent   ! note: this must include the faces and areas of Cartesian cells \\
%\indent   TYPE(CUTCELL\_LINKED\_LIST\_TYPE), POINTER :: CUTCELL\_FACE\_LIST \\
%\noindent END TYPE CUTCELL\_TYPE \\
%
%\noindent TYPE FACET\_TYPE \\
%\indent   INTEGER :: VERTEX(3)=0,SURF\_INDEX=0,BOUNDARY\_TYPE \\
%\indent   CHARACTER(LABEL\_LENGTH) :: SURF\_ID='null' \\
%\indent   REAL(EB) :: NVEC(3)=0.\_EB,AW,EW,KW,DN,RDN,RHO\_F,U\_TAU,Y\_PLUS,TMP\_F, \\
%\indent TMP\_G,QCONF,HEAT\_TRANS\_COEF,QRADIN,QRADOUT \\
%\indent   REAL(EB), ALLOCATABLE, DIMENSION(:) :: RHODW,ZZ\_F \\
%\indent   REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: ILW \\
%\indent   TYPE(CUTCELL\_LINKED\_LIST\_TYPE), POINTER :: CUTCELL\_LIST \\
%\noindent END TYPE FACET\_TYPE \\
%
%\noindent TYPE(FACET\_TYPE), ALLOCATABLE, TARGET, DIMENSION(:) :: FACET \\
%
%\noindent TYPE VOLUME\_TYPE \\
%\indent   INTEGER :: VERTEX(4)=0 \\
%\indent   CHARACTER(LABEL\_LENGTH) :: MATL\_ID='null' \\
%\noindent END TYPE VOLUME\_TYPE \\
%
%\noindent TYPE(VOLUME\_TYPE), ALLOCATABLE, TARGET, DIMENSION(:) :: VOLUME \\
%
%\end{myfont}
%
%
%\subsubsection*{Proposed Data structures for the FDS CCIBM code:}
%
%The previously defined derived type \texttt{CUTCELL\_LINKED\_LIST\_TYPE} refers to \texttt{INBOUNDARY} cut-faces as defined in the Matlab data structures. Also, previously there was no need for a proper volume cut-cell linked list type, because (from what I understand) things on previous development didn't get to this stage (looks like Charles Luo was playing with one CC on the INIT\_IBM subroutine, line $\simeq 2750$ of \texttt{geom.f90}).
%We modify the previously derived types keeping the data already defined and needed, and considering the Matlab code data containers:
%\newline
%
%\noindent  - One type \texttt{CUTFACE\_LINKED\_LIST\_TYPE} to create a linked list with the data corresponding to the \texttt{CUT\_FACE} data structure in the Matlab code. The type proposed is:
%
%\begin{myfont}
%\noindent TYPE CUTFACE\_LINKED\_LIST\_TYPE \\
%\indent  INTEGER :: INDEX                                                                        ! Cut- face index \\
%\indent  TYPE(CUTFACE\_LINKED\_LIST\_TYPE), POINTER :: NEXT    ! nxt cut-face \\
%\indent   REAL(EB), DIMENSION(NDIM) :: XYZCEN                                 ! cut-face centroid coords \\
%\indent   INTEGER :: NVERTEX, STATUS, IJK, AXIS, BCTYPE(NSCALARS)            \\
%\indent   ! Vertex coordinates, to be allocated (NDIM,NVERTEX) \\
%\indent   REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: XYZVERT \\
%\indent   REAL(EB) :: BCVAL(NSCALARS) \\
%\noindent END TYPE CUTFACE\_LINKED\_LIST\_TYPE \\
%\end{myfont}
%Other fields can be added as needed.
%
%\noindent  - One new type \texttt{CUTCELL\_LINKED\_LIST\_TYPE} to create a linked list with the data corresponding to the \texttt{CUT\_CELL} data with basic cut-cell info for flux definitions and volume integrals:
%\begin{myfont}
%\noindent TYPE CUTCELL\_LINKED\_LIST\_TYPE \\
%\indent  INTEGER :: INDEX                                                                        ! Cut-cell index \\
%\indent  TYPE(CUTCELL\_LINKED\_LIST\_TYPE), POINTER :: NEXT    ! nxt cut-cell \\
%\indent  REAL(EB) :: VOL                                                                           ! cut-cell volume \\
%\indent  REAL(EB), DIMENSION(NDIM) :: XYZCEN                                 ! cut-cellcentroid coords \\
%\indent  ! List of cut-faces that are boundary of this cut-cell \\
%\indent  TYPE(CUTFACE\_LINKED\_LIST\_TYPE), POINTER :: CUTCELL\_CUTFACES\_LIST \\
%\indent  INTEGER :: IJK, HINTFLG ! IJK underlying Cartesian cell. \\
%\noindent END TYPE CUTFACE\_LINKED\_LIST\_TYPE \\
%\end{myfont}
%Other fields added as needed (\texttt{RHO,TMP,DIV,ZZ} etc.).
%
%\noindent We need to change the name of the linked list to cut-faces in the \texttt{FACET\_TYPE} definition: \\
%\begin{myfont}
%\indent   TYPE(CUTCELL\_LINKED\_LIST\_TYPE), POINTER :: CUTCELL\_LIST \\
%\end{myfont}
%\noindent  to
%\begin{myfont}
%TYPE(CUTFACE\_LINKED\_LIST\_TYPE), POINTER :: FACET\_CUTFACES\_LIST \\
%\end{myfont}
%\noindent to define correctly the type of geometry entity the list refers to.
%
%Then we would define: \\
%\begin{myfont}
%\noindent !LL of cut-faces for mesh M. \\
%\noindent TYPE(CUTFACE\_LINKED\_LIST\_TYPE), POINTER :: CUT\_FACE\_LL \\
%\noindent !LL of cut-cells for mesh M. \\
%\noindent TYPE(CUTCELL\_LINKED\_LIST\_TYPE), POINTER :: CUT\_CELL\_LL \\
%\end{myfont}
%
%From my understanding on linked lists, using a \texttt{CUTCELL\_CUTFACES\_LIST} for each cut-cell defined in \texttt{CUT\_FACE\_LL, CUT\_CELL\_LL} and also for some \texttt{FACET} would mean we have to replicate cut-face data at least 3 times? Cut-faces stored on the different lists would be individualized by their global cut-face \texttt{INDEX}. This seems convoluted and prone to bugs.
%The other option is to use the \texttt{CUT\_FACE\_LL} linked list to generate an allocatable data structure \texttt{CUT\_FACE} with same fields, and define also \texttt{CUT\_CELL} in the same manner. The \texttt{CUT\_CELL} data structure size can be estimated from the beginning knowing the number of Cartesian cells that are cut. Then, \texttt{FACET} would have an allocatable index array to the \texttt{CUT\_FACE} data structure.
%
%\textit{We need to discuss this..}
%
%
%\subsubsection*{Containers for Pressure and scalar unknowns numbering:}
%
%\textit{To come..}
%
%\newpage
%
%
%
%
%\begin{figure}
%\begin{center}
%\begin{tikzpicture}[scale=1]
%
%\coordinate (O) at (0,0);
%\coordinate (P) at (10,10);
%
%\draw[help lines] (O) grid [step=2cm] (P);
%
%\end{tikzpicture}
%\caption{Structured, staggered grid in FDS.  Mesh variables belong to {\ct M=>MESHES(NM)}. }
%\label{fig:struc_stag_grid}
%\end{center}
%\end{figure}
%
%
%\begin{figure}
%\begin{center}
%
%\newcommand{\Depth}{3}
%\newcommand{\Height}{3}
%\newcommand{\Width}{3}
%
%\tdplotsetmaincoords{70}{30}
%
%\begin{tikzpicture}[scale=2,tdplot_main_coords]
%
%  \coordinate (O) at (0,0,0);
%  \coordinate (A) at (0,\Width,0);
%  \coordinate (B) at (0,\Width,\Height);
%  \coordinate (C) at (0,0,\Height);
%  \coordinate (D) at (\Depth,0,0);
%  \coordinate (E) at (\Depth,\Width,0);
%  \coordinate (F) at (\Depth,\Width,\Height);
%  \coordinate (G) at (\Depth,0,\Height);
%
%  \tdplotsetcoord{P}{1}{45}{45}
%  \draw[ultra thick,->] (0,0,0) -- (1,0,0) node[anchor=north east]{$x$};
%  \draw[ultra thick,->] (0,0,0) -- (0,1,0) node[anchor=north west]{$y$};
%  \draw[ultra thick,->] (0,0,0) -- (0,0,1) node[anchor=north east]{$z$};
%%  \draw[-stealth,color=red] (O) -- (P);
%%  \draw[dashed, color=red] (O) -- (Px);
%%  \draw[dashed, color=red] (O) -- (Py);
%%  \draw[dashed, color=red] (O) -- (Pz);
%%  \draw[dashed, color=red] (Px) -- (Pxy);
%%  \draw[dashed, color=red] (Py) -- (Pxy);
%%  \draw[dashed, color=red] (Px) -- (Pxz);
%%  \draw[dashed, color=red] (Pz) -- (Pxz);
%%  \draw[dashed, color=red] (Py) -- (Pyz);
%%  \draw[dashed, color=red] (Pz) -- (Pyz);
%%  \draw[dashed, color=red] (Pxy) -- (P);
%%  \draw[dashed, color=red] (Pxz) -- (P);
%%  \draw[dashed, color=red] (Pyz) -- (P);
%
%  \draw[fill=gray!50,opacity=0.5] (O) -- (A) -- (E) -- (D) -- cycle;
%  \draw[fill=gray!50,opacity=0.5] (O) -- (A) -- (B) -- (C) -- cycle;
%  \draw[fill=gray!50,opacity=0.5] (A) -- (B) -- (F) -- (E) -- cycle;
%
%  \draw[fill=gray!50,opacity=0.5,thick] (D) -- (E) -- (F) -- (G) -- cycle;
%  \draw[fill=gray!50,opacity=0.5,thick] (C) -- (B) -- (F) -- (G) -- cycle;
%  \draw[fill=gray!50,opacity=0.5,thick] (O) -- (C) -- (G) -- (D) -- cycle;
%
%  %% Following is for debugging purposes so you can see where the points are
%  % These are last so that they show up on top
%  \foreach \xy in {O, A, B, C, D, E, F, G}{
%      \node at (\xy) {\xy};
%  }
%
%  % DX(I)
%  \draw [blue,<-] (0,-.4*\Width,0) -- (.35*\Depth,-.4*\Width,0);
%  \draw [blue,->] (.65*\Depth,-.4*\Width,0) -- (\Depth,-.4*\Width,0);
%  \draw [blue,-] (0,-.5*\Width,0) -- (0,-.1*\Width,0);
%  \draw [blue,-] (\Depth,-.5*\Width,0) -- (\Depth,-.1*\Width,0);
%  \draw (.5*\Depth,-.4*\Width,0) node {\ct DX(I)};
%
%  % DY(J)
%  \draw [blue,<-] (1.3*\Depth,0,0) -- (1.3*\Depth,0.35*\Width,0);
%  \draw [blue,->] (1.3*\Depth,0.65*\Width,0) -- (1.3*\Depth,\Width,0);
%  \draw [blue,-] (1.1*\Depth,0,0) -- (1.4*\Depth,0,0);
%  \draw [blue,-] (1.1*\Depth,\Width,0) -- (1.4*\Depth,\Width,0);
%  \draw (1.3*\Depth,.5*\Width,0) node {\ct DY(J)};
%
%  % DZ(K)
%  \draw [blue,<-] (1.3*\Depth,\Width,0) -- (1.3*\Depth,\Width,0.45*\Height);
%  \draw [blue,->] (1.3*\Depth,\Width,0.55*\Height) -- (1.3*\Depth,\Width,\Height);
%  \draw [blue,-] (1.1*\Depth,\Width,\Height) -- (1.4*\Depth,\Width,\Height);
%  \draw (1.3*\Depth,\Width,.5*\Height) node {\ct DZ(K)};
%
%  % XC(I),YC(J),ZC(K)
%  \filldraw [black] (.5*\Depth,.5*\Width,.5*\Height) circle (2pt);
%  \node at (.5*\Depth,1.3*\Width,1.2*\Height) {\ct XC(I),YC(J),ZC(K)};
%  \node at (.5*\Depth,1.3*\Width,1.3*\Height) {\ct RHO(I,J,K),ZZ(I,J,K,N),TMP(I,J,K),H(I,J,K)};
%  \draw [blue] (.5*\Depth,.55*\Width,.55*\Height) -- (.5*\Depth,1.25*\Width,1.15*\Height);
%
%  % velocity components
%  \draw [->, ultra thick, >=latex] (.9*\Depth,.5*\Width,.5*\Height) -- (1.1*\Depth,.5*\Width,.5*\Height);
%  \draw [->, ultra thick, >=latex] (.5*\Depth,.9*\Width,.5*\Height) -- (.5*\Depth,1.1*\Width,.5*\Height);
%  \draw [->, ultra thick, >=latex] (.5*\Depth,.5*\Width,.9*\Height) -- (.5*\Depth,.5*\Width,1.1*\Height);
%
%  \node at (\Depth,.5*\Width,.4*\Height) {\ct U(I,J,K)};
%  \node at (.75*\Depth,\Width,.5*\Height) {\ct V(I,J,K)};
%  \node at (.45*\Depth,.6*\Width,1.1*\Height) {\ct W(I,J,K)};
%
%  % shade solid volume
%  \coordinate (V1) at (0,.75*\Width,\Height);
%  \coordinate (V2) at (.75*\Depth,0,\Height);
%  \coordinate (V3) at (0,0,.1*\Height);
%
%  % face normal
%  \coordinate (N1) at (.25*\Depth,.25*\Width,.7*\Height);
%  \coordinate (N2) at (.4*\Depth,.4*\Width,.55*\Height);
%  \draw [->, thick, >=latex] (N1) -- (N2);
%  \draw [fill=red,opacity=0.5,thick] (V1) -- (V2) -- (V3) -- cycle;
%  \draw [fill=green,opacity=0.5] (V2) -- (V1) -- (C) -- (V3) -- cycle;
%
%  % facet vertices intersecting with box
%  \draw (C) -- (G);
%  \filldraw [gray] (V1) circle (2pt);
%  \filldraw [gray] (V2) circle (2pt);
%  \filldraw [gray] (V3) circle (2pt);
%  \node at (.75*\Depth,0,.9*\Height) {$V_2$};
%  \node at (-.05*\Depth,.75*\Width,1.05*\Height) {$V_1$};
%  \node at (-.05*\Depth,-.05*\Width,.1*\Height) {$V_3$};
%
%\end{tikzpicture}
%
%%\draw (0.30,0.25) node {$\mathbf{v}_1$};
%%\draw (4.80,3.10) node {$\mathbf{v}_2$};
%%\draw (2.20,4.00) node {$\mathbf{v}_3$};
%%
%%\draw [->, very thick, >=latex] (-2.5,5.0) -- (-1.5,5.0);
%%\draw [->, very thick, >=latex] (-2.5,3.0) -- (-1.5,3.0);
%%\draw [->, very thick, >=latex] (-0.5,5.0) -- (0.5,5.0);
%%\draw [->, very thick, >=latex] (-0.5,3.0) -- (0.5,3.0);
%%\draw (0.5,3.3) node {$\mathbf{x}_{velo}$};
%%
%%\draw [->, very thick, >=latex] (0.2,0.5) -- (0.0,3.0);
%%\draw (-0.15,1.8) node {$\mathbf{r}$};
%%
%%\draw [->, very thick, >=latex] (0.2,0.5) -- (-1,2.0);
%%\draw (-0.5,1.0) node {$\mathbf{n}$};
%%
%%\draw [->, >=*] (1.1,1.545) -- (-1.3,4.7);
%%\draw (-1.5,4.3) node {$\mathbf{x}_{int}$};
%%
%%\draw [|-] (0.9,1.37) -- (0.57,1.8);
%%\draw [-|] (0.24,2.25) -- (-0.19,2.85);
%%\draw (.4,2.02) node {$\delta n$};
%%
%%\draw [->, very thick, >=latex] (1.08,1.55) -- (2.35,1.55);
%%\draw [->, very thick, >=latex] (1.1,1.545) -- (0.45,2.4);
%%\draw [->, very thick, >=latex] (1.1,1.545) -- (1.75,2.6);
%%\draw (2.9,1.55) node {$\mathbf{\hat{s}} \, (\mathbf{\bar{e}}_1)$};
%%\draw (1.9,2.60) node {$\mathbf{\hat{t}}$};
%%\draw (0.6,2.60) node {$\mathbf{\hat{n}}$};
%%
%%\draw [->, very thick, >=latex] (-2,0) -- (-2,1);
%%\draw [->, very thick, >=latex] (-2,0) -- (-1,0);
%%\draw (-0.9,-0.25) node {$\mathbf{e}_1$};
%%\draw (-2.35,1) node {$\mathbf{e}_3$};
%%
%%\draw (1,4.7) .. controls (2.5,2) .. (1.1,1.545);
%%\draw [->] (-1.15,4.500) -- (1.10,4.500);
%%\draw [->] (-0.58,3.750) -- (1.52,3.750);
%%\draw [->] ( 0.20,3.003) -- (1.90,3.003);
%%\draw [->] ( 0.54,2.280) -- (2.15,2.280);
%%\draw [->] ( 0.90,1.800) -- (1.80,1.800);
%%\draw (1.5,5.0) node {$f(n)$};
%
%\caption{Cartesian cell {\ct (I,J,K)} in FDS.}
%\label{fig:cartesian_cell}
%\end{center}
%\end{figure}

\bibliography{../Bibliography/FDS_general,../Bibliography/FDS_refs,../Bibliography/FDS_mathcomp,../Bibliography/sv_fire,../Bibliography/sv_graphics}

\addcontentsline{toc}{chapter}{References}

\end{document}
